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1. INTRODUCTION 

Future research development, as well as production activities in space by 
U.S. federal and private agencies, will depend on the Space Shuttle and its 
derivative versions as a principal space transportation system. This depend- 
ence requires improved designs or techniques to extend the life, upgrade 
performance, reduce weight, lower operational costs, and generally improve the 
functional capability of the main propulsion system. The engines for this main 
propulsion system are advanced high pressure engines operating on oxygen and 
hydrogen. A need therefore exists to investigate, develop, and define basic 
concepts in support of the main propulsion system improvements. One basic 
area that has been investigated is the hot gas flow nonuniformities that occur 
within the manifold, duct work, and main injector. Nonuniformities result from 
highly distorted and mismatched flows within the ducts which create severe 
environments for the system components, thus limiting their useful life. 

Development and verification tests of the Space Shuttle Main Engine (SSME) 
and results of computational modeling have shown that the three gas transfer 
tubes have an uneven flow distribution with large areas of separated flow. The 
outer transfer tubes each carry approximately twice the amount of gas as the 
center tube. This causes the energy of the gas to be much higher in the outer 
tubes. Flow from the tubes impinge upon the main injector liquid oxygen posts 
which bend under the static load of the gas flow, the bending being more pro- 
nounced in line with the outer tubes. To alleviate this phenomenon and to 
keep the posts cooler, shields, linking pairs of posts in the outermost row, 
were installed. The design alteration enhanced the injector life, however, 

LOX post failures show that this change alone is insufficient to grant speci- 
fied life at equal to or greater rated power levels. Incorporation of the 
shields also affects circumferential flow in the annulus and degrades the 
engine performance and necessitates higher operating temperatures in the 
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turbines. To improve the SSNE design and for future use in the design of new 
generation rocket engines, more experimental data are needed. These data, 
combined with previous test data and integrated into analytical methods, can 
be used to establish a base for design of high energy hot gas flow systems. 

In addition the prebumers of the SSME still experience serious problems 
related also to the cracking of LOX posts and the frisbee in the inlet mani- 
folds. Operation of the fuel prebumer creates flows which cause turbine 
blade cracking in the high pressure fuel turbopump. The most likely cause of 
these problems is the extreme thermal environments caused by start-up and 
shut-down of the engine. A detailed transient analysis of these flows is 
necessary in order to help alleviate these problems. 

This report describes results of efforts by personnel of the Computational 
Mechanics Group at the Lockheed-Huntsville Engineering Center to assist the 
computational staff of NASA-MSFC in developing analytical capabilities for 
modeling hot gas flow on the fuel side of the SSME. The primary objective was 
to develop and deliver a computer code Which produces a computational grid for 
the three-duct Hot Gas Manifold (HGM) on the fuel side of the SSME. A 
secondary objective was to provide input for an accurate computation of the 

v, 

flow characteristics inside the fuel side prebumer during the start-up 
transient of the SSME. 
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2. TECHNICAL APPROACH 
2.1 HGM COMPUTATIONAL GRID CODE 


Nearly all computational codes which are available to numerically solve 
the three-dimensional fluid flow equations are designed to be applied to a 
well structured grid model of the flow region. To perform the calculations, 
generalized independent variables are introduced which transform the physical 
coordinates, (x, y, z), into general curvilinear coordinates, (n_, n 3 )- 

Thus the physical domain must be gridded as a single or series of hexahedral 
zones described by eight comer points, 12 edges, and six surfaces. Such an 
arbitrary zone is shown in Fig. 2-1. 



Fig. 2-1 Hexahedral Element Showing Local Intrinsic Coordinates 
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The approach used in the current study was to provide an algebraic grid 
generation code which would produce the Cartesian coordinates (x, y, z), for 
points along the lines of constant n^, and Basic mathematical tech- 

niques taken from analytic geometry and vector algebra were employed to des- 
cribe a hexahedral zone in terms of piecewise continuous analytic functions 
which represent the zonal edges and surfaces. 


An intrinsic curvilinear coordinate system can be produced by mapping a 
unit cube onto the simply connected hexahedral zone. What is needed is a 
transformation function that will map a unit cube in (n. , h,) space uni- 

valently onto the hexahedral volume of interest thus producing the required 
intrinsic coordinate system. A procedure which produces the desired result is 
referred to as either the method of transfinite interpolation or multi-variate 
blending function, interpolation (Ref. 1). A brief description of this method 
follows . 


Let F(n ; n 2> t^) be a vector-valued functional representing the region R 
of interest in curvilinear space. Then as n 3 > range over R, F traces 

out the region in Euclidean space (x, y, z). Also let 4>, and V be 
blending functions which obey the cardinality conditions: 

4> i (h 1 =l) 

(T» 2 =m) 

Vv n) 


| 1 if i = 1 
(o if i # 1 

1 1 if j = m 
0 if j / m 
j 1 if k = n 
I 0 if k y n 


The simplest form of blending functions meeting these conditions are 


(n. ) = l - n 

O 1 l 

W = 1 - ^2 

X (Hi = 1 - 


= n 1 

VV ’ Y 
W ' " 3 


Then a trilinearly blended interpolant of F, which will map a unit cube onto 
the region R is given by 
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U(n i' n 2' 


where 




■x(n 1 ,n 2 ,n ;J )‘ 

y ( WV 

. Z (n :L .Ti 2 ,n 3 )J 


= (i-n 1 )F(o n 2 ,n 3 ) + ^ f<i . n 2 ,n 3 > 

+ (i-n 2 )F(n lf o,n 3 > + t^fo^.i.Hj) 


+ (l-n 3 )F(n lf n 2 ,0) + T1 3 F(T1 1 » T1 2 ,1) 

- (l-n )(i-n 2 >)F(0,o,n ) - (i-n 1 )n 2 F(o, i.n 3 > 


- n 1 (l-n 2 )F(i,o,n 3 ) - ly^FU.l.iij) 

- d-n 1 )n 3 F(0,n 2 ,0) - (i-n 1 )n ;} F(0,n 2 ,i) 

- n 1 (i-n 3 )F(i,n 2 ,0) - n^FU.iyl) 

- (l-n 2 )(l-n 3 )F(n 1 ,0,0) - (i-n 2 )n 3 F(n 1 o,i) 

- n 2 (i-n ;j )F(n 1 ,l,0) - T ’ 2 I1 3 F{T, i ,1,1) 

+ (i-n 1 )(i-n 2 )(i-n 3 )F(0,0 l 0) + (i-n^(l-n 2 )n 3 F(0,0,i) 
+ . (i-h 1 )n 2 (i-n 3 )F(0,i,o) + (i-n 1 )ti 2 n 3 F(0,i,D 
+ n 1 (i-n 2 )(i-n 3 )F(i t 0,0) + n 1 (i-n 2 >n 3 F(l,0,l) 

+ T1 1 n 2 <1_n 3 ) F ( 1,1,0 > + T, 1 I, 2 n 3 F(1,1,1) 


F(0,n 2 ,n 3 ) represents a surface with n 3 =0, etc. 
F(0,0,n 3 > represents an edge with n =n 2 =0, etc. 
F(0 ,0,0) represents a point with n 3 =n 2 =n 3 =0, etc. 
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In two dimensions this equation performs a bilinear interpolation over an 
arbitrary region consisting of four distinct comers simply connected by four 
edges, where 


F(0,n 2 > = 


along EDGE. 

4 

F ( 0 , 0 ) 


at 

POINT. 

m,n 2 ) = 


along EDGE 2 

F(0,1) 


at 

POINT 

F(n 1 ,0) = 

*;> 

along EDGE^ 

F(1,0) 

- 

at 

POINT 

F(n 1 ,D = 


along EDGE 3 

F(l,l) 

-S' 

at 

POINT 


and the interpolation equation for F could be rewritten as 

U = [*] = (1-H 1 )EDGE 4 + n 1 EDGE 2 + (l-n 2 >EDGE 1 + n 2 EDGE 3 

- (i-n 1 )(i-n 2 )POiMT 1 - (l-n 1 )n 2 POiNT 4 

- n 1 (i-n 2 )POiNT 2 - n 1 n 2 POiNT 3 

In this equation EDGE, represents a vector-valued functional along edge 

4 

four, etc. Examination of this equation shows that it performs linear 
interpolations between EDGE^ and EDGE 3 and between EDGE 2 and EDGE 4> hence the 
term bilinear interpolation. Hence, if the functional for each edge can be 
derived and is analytic, a grid or mesh of node points can be generated by 
substituting values of and ti 2 into the equation. 

In three dimensions the general equation above performs trilinearly 
blended interpolation over an arbitrary region consisting of eight distinct 
corner points simply connected by 12 edges, where 
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F(0,n 2 ,n 3 ) 

F(0,o,n 3 ) 

F(0 ,0,0) 

and which can be rewritten as 
U = (jjJ 

. (i-n . )side + n,siDE. t (l-n.)SiOE, + h.side. + (i-h,)side, + n side ' 

1516 2 224 3 133 

- (l-n.. )(l-n,)EDGE - (l-n.)n-EDGE 0 - n. (l-n.)EDGE. - n h-EDGE 

1 2 5 1 2 8 1 2 6 1 2 7 

- u-n 1 )(l-n 3 )EDGE 4 - (i-n 1 )Ti 3 EDGE 12 - n 1 (i-n 3 )EDGE 2 - dredge 

- (l-n i )Cl-n 3 )EDGE 1 - (l-n 2 )n 3 EDGE g - Tl 2 (l-n 3 )EDGE 3 - T ’ 2 n 3 EDGE ll 

+ (l-n-)(l-O<l-OP0Iirr + (l-iv)(l-n)n POINT. + (l~n )T»„(i-n,) point 
1 2 3 1 1 23 5 123 4 

+ (i-n 1 )n 2 n 3 POiNT 8 + n 1 (l-n 2 )(i-n 3 )POiNT 2 + r^U-n^t^POiNT^ 

+ rt l n 2 (i-n 3 )POiNT 3 + r, 1 n 2 n 3 POINT 7 
where SIDE., EDGE., POINT represent vector-valued functionals 

13 K 

on the surfaces, along the edges, and at the corner points, 
respectively. 
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This equation reduces to the previous two-dimensional analog along any 
flat surface or along any surface in which a straight line can be drawn 
between any two opposing edges such that the line lies entirely within the 
surface. 

With the general transformation, any point in local coordinates n 2 » n 3 
can be related to the physical Cartesian coordinates x, y, z. The entire grid 
of discrete points is generated in the HGM code using this concept. This 
general interpolant can accommodate any stretching function for concentrating 
points near walls or regions of large gradients. Furthermore, the edges of 
the hexahedral can be segmented allowing another means of grid spacing control. 

2.2 SSME FUEL-SIDE PREBURNER ANALYSIS 


Operation of the fuel prebumer is initiated by opening propellant 
control yalves and igniting the flow with the augmented spark ignited (ASI). 
The configuration of the FPB is shown in Fig. 2-2. The flow passage from the 
fuel prebumer oxidizer valve (FPOV) into the oxidizer manifold above the 
interpropellant plate is shown in this figure. The oxidizer flow to the ASI 
is supplied from a bleed line from the FPOV. The valves which control this 
flow are shown schematically in Fig. 2-3; note that there is not a separate 
fuel control valve for each prebumer. Flows from the fuel and oxidizer 
manifolds enter the FPB combustor through the injector face plate, one third 
of the layout of which is shown in Fig. 2-2. The prebumer liner extension 
shown in Fig. 2-2 slides inside of the sealer groove as shown in the upper 
left and right sides of Fig. 2-5. Thus, the prebumed gases flow, via the 
dome at the bottom of the chamber, onto the turbine blades, then following the 
path through the 180-deg turnaround duct to the fuel bowl and transfer ducts 
of the HGM and onto the LOX posts and into the main combustion chamber. 

The presence of the baffles shown in Figs. 2-2 and 2-4 create symmetry 
planes for computational analysis but remove the simplicity of axis symmetry. 
In addition, any model geometry should include the influence of the large nut 
at the bottom of the prebumer chamber which holds the dome in place. 
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Fig. 2-2 Fuel Prebumer Configuration 
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Fig. 2-4 Face Plate Showing One-Third of Preburner 

Distribution of Injection Holes and Elements 
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Fig. 2-5 Cross-Section of Lower Half Fuel Side Prebumer Chamber 
Indicating at Top Left and Right Procedure for Insertion 
of Upper Preburner Liner Extension 
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Our approach is to use the symmetry created by the baffles and model only 
one-third of the geometry, and the large nut on top of the dome should be in- 
corporated into the computational grid. Furthermore, in order that the most 
accurate analysis of the prebumer be made for portions of the startup 
transient stage of its operation, a thorough study was made of the average 
chamber parameters for the first few seconds of operation. Data were obtained 
from the latest Digital Transient Model (DTM) results. The DTM is a one- 
dimensional lumped parameter model of the entire SSME which was originally 
developed by Rocketdyne and modified by D.C. Seymour of the Performance 
Analysis Branch of the Propulsion Lab at NASA-MSFC, to more closely match 
actual engine test data obtained at MSFC. These data were used to obtain, as 
accurate as possible, the face plate inlet conditions and downstream chamber 
exit conditions for use in the boundary treatment of a computational analysis. 
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3 . RESULTS 

3.1 HGM COMPUTATIONAL GRID CODE 

3.1.1 HGM Geometry 

A schematic representation showing the construction of the geometry which 
must be modeled is presented in Fig. 3-1. The 180-deg turnaround duct is dis- 
played in a cutaway fashion to show the turn. Only half of the manifold is 
represented since, due to the plane of symmetry which divides the center 
transfer duct in half, only half need be computationally modeled. 



Fig. 3-1 Schematic Showing Components of HGM Geometry to be Modeled 
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The geometry was generated in five separate pieces or zones. The five 
zones are shown in Fig. 3-2. Each zone has a general hexahedral shape and in 
Fig. 3-2 most of the eight comer points of each are clearly marked. Figure 
3-3 shows the Cartesian coordinate system relative to which the position of 
each node in the grid is referenced. Except for the hole perimeters in the 
outer wall of Zone 3 and the fairing at the entrances to Zones 4 and 5 the 
edges of each zone can be described as piecewise continuous segments composed 
of either straight lines or circular arcs. In addition, excluding the two 
special regions previously mentioned, the surfaces of each zone can be gen- 
erated by rotating an edge about the appropriate axis. For Zones 1, 2, and 3 
the X axis is the axis of revolution. 

3.1.2 HGM Grid 

Zone 1, the turnaround duct (TAD), is composed of 94 nodes in the stream- 
wise direction, 72 nodes in the circumferential direction (0 to 180 deg), and 
21 nodes across the duct between inner and outer surfaces. Two perspectives 
showing this nodal distribution are shown in Fig. 3-4. 

Zones 2 and 3 make up the bowl section of the HGM and contain 59 nodes 
in the x-direction, from bowl entrance to rear of the bowl, 107 nodes in the 
circumferential direction and 21 nodes between the inner and outer surfaces. 
Distribution of nodes in the bowl is presented in Figs. 3-5 through 3-7. 

Zones 4 and 5 comprise the right and half of the middle transfer duct 
portions of the HGM. The right duct has been generated with 29 nodes along 
the duct axis and 27 x 20 nodes in a cross section. The similar distribution 
in Zone 5 is 29 x 27 x 13. Surface and cross-section grids for these two 
zones are given in Figs. 3-8 and 3-9. A closeup of the surface mesh in the 
vicinity of the fairing for each transfer duct is contained in Fig. 3-10. 
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Fig. 3-4 Grid Plots of TAD Internal Grid (Bottom) and of Part 
of the Grid Distribution on Inner and Outer Surfaces 
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Fig. 3-5 Grid Plots for Inner and Outer Surfaces of Zones 2 and 3 

as Well as Internal Grid at the Common Plane of Intersection 
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Fig. 3-6 Grid Plot Showing Node Distribution in Zones 2, 3, 4, and 5 on 
a Plane Which Passes Through the Center of Both Transfer Ducts 
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Fig. 3-7 Grid Plot of Computational Mesh on Outer Surface of Zone 3 
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Fig. 3-9 Middle Transfer Duct Surface 

and Cross-Section Nodal Distribution 
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Fig. 3-10 Expanded View of Surface Mesh for Region 
of Fairing on Both Transfer Ducts 
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The internal grid resolution presented in the previous figures is ^ 
adequate for a laminar viscous computation but would need to be modified for 
application to a turbulent computation. The procedure for doing this will be 
described in Section 3.1.4. 

3.1.3 Computer Code 


A source listing of the three-duct HGM geometry code is provided in 
Appendix C. A concise input guide and an input listing for the grids 
displayed in Figs. 3-4 through 3-9 is also included in Appendixes A and B, 
respectively. 


In Fig. 3-11 a primary calling sequence flow chart is shown. A brief 
explanation of the function of each subroutine follows: 


• INITIAL - Reads the first two lines of the input file, initializes 

coefficient arrays, and defines logical unit numbers and 
counters . 

• INPUT - Reads the remainder of the input file and sets all parameters 

to be used in remaining subroutines. 

• GRID - The controlling subroutine for the generation of the spatial 

coordinates of each node in each zone. 

• ETABC - Calculates the values of tjj_, Ti2* TI3 along the I, J, and K 

directions for each hexahedral shaped section of each zone. 

• EDGE - Determines the Cartesian coordinates for nodes along each edge 

of each side of each section of each zone using the bilinear/ 
trilinear interpolation scheme. 

• SURFACE - Determines the Cartesian coordinates of nodes on a three- 

dimensional surface using the trilinear interpolation scheme. 
Here, all outer and internal surface nodes are calculated 
from the previously determined edge distribution. 

• OUTPUT - Provides printed output and stores geometry in File 20 for 

use in plotting or as input to an integration code. 
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The output to File 20 is in the format to be input as a multi-grid 
geometry file to the PL0T3D plotting code. 

The code is written in standard FORTRAN. However, the listing provided 
in Appendix C is a VAX 11/785 version and could contain some generic VAX 
statements which would have to be translated if used on other than a VAX. 

3.1.4 Code Implementation 

To implement the HGM geometry code as it appears in the listing of 
Appendix C, the input listing of Appendix B must correspond to logical Unit 5. 
Unit 6 must be assigned to the written output and Unit 20 to the geometry file 
which will contain the x, y, z coordinates for each node in each zone. 

Each separate input file is labeled to indicate the zone being described 
by the following card images. These zone labels correspond to those shown in 
Fig. 3-2. Because of the number of nodes in each zone it is recommended that 
each zone be run separately. In fact as shown in Appendix B, each zone input 
is a separate file since all contain card types 1, 2, and 3. If more than one 
zone is to be run together (e.g., Zones 2 and 3), only one set of card types 
1, 2, and 3 should appear at the top of the input file. 

A detailed description of the input file to the code is provided in 
Appendix A. Modifications to the geometry can be facilitated by studying the 
input guide while observing both the input files and the detailed grid pictures 
presented in Figs. 3-3 through 3-10. Redistributions of nodes can be accomp- 
lished by making minor modifications to the input files. For example for a 
turbulent computation if the nodes in Zone 4 near the wall require redistribu- 
tion closer to the wall then a change would need to be made to card type 9 on 
line 11 of the input file for that zone. The 7.0 appearing in the and n 3 
positions could be changed to 10.0 (see page A-10). If more nodes were 
desired in a cross plane then the 27 and 29 on card type 8 could both be 
increased. 
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Dimensioning in the program has been kept to a minimum. The largest 
dimensioned arrays in the bulk of the code are NOD ENUM( 10000) and 
N0DE(5 ,10000) . At the very end of all computations the PL0T3D file is gen- 
erated. In this subroutine the x, y, z coordinate arrays are each dimensioned 
to 155,000. The 10,000 corresponds to twice the maximum number of nodes in a 
plane perpendicular to the marching direction (n direction input on card 7) 
for creating the geometry. The 155,000 must be equal to or greater than the 
number of nodes in the largest zone, which is the TAD. In Zones 1, 2, and 3 
the direction of is from TAD entrance to bowl back wall; the direction 
is from inner to outer wall; and is directed from side opposite transfer 
ducts circumferentially. In zones 4 and 5, increases in the streamwise 
direction from bowl outward, and x form the cross planes in each duct. 

Note that the code is designed to output each zone of the HGM so that 
each has one cross plane in common with the preceding zone(s). This must be 
remembered for incorporating the grid into a flowfield solver code. The 
geometry must be integrated in a multi-block or multi-zone fashion. If the 
computer available has large enough core memory or if it is a large virtual 
machine then Zones 2 and 3 can easily be combined into one larger zone since 
at all common planes the nodal positions match exactly. However, the TAD zone 
matches exactly to Zone 2 but only every other circumferential plane in Zone 3 
matches with a corresponding plane in the TAD. This mismatch can be corrected 
by making appropriate minor modifications to the input file for Zone 1. 


3.2 PREBURNER ANALYSIS 


An accurate analysis of the fluid dynamic environment in the prebumer 
during the startup transient portion of its operation demands careful treat- 
ment of the inflow and outflow boundary conditions as well as the most up-to- 
date chemistry model for the reacting components in the flow. Results of a 
comprehensive study of these two items is presented in the following two 
sections . 
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3.2.1 Inflow and Outflow Boundary Treatment 

Data were obtained in the form of printed output from the most recent DTK 
analysis of the SSHE startup transient conditions. This included injector fuel 
and oxidizer flow rates, ignition fuel, and oxidizer flow rates. From this, 
average chamber temperature and static pressure data tables were generated. 
These data are displayed graphically in Figs. 3-12 through 3-16. At any 
instant in time during the 5-sec startup process, these data provide informa- 
tion necessary to specify mass and composition flow rates at the inlet plane 
as well as downstream static pressure at the exit plane for performing a flow- 
field computation. 

During any selected time interval, downstream static pressure would be 
obtained from a table look-up of the data shown in Fig. 3-16. Total flow 
rates are obtainable from the remaining data, but the specific strategy for 
distributing the mass flow over the face plate inlet plane requires careful 
analyses. Observation of Fig. 2-4 reveals that there is clearly a nonuniform 
distribution of mass injection. A practical approach which allows for a 
reasonable approximation to the physical situation would be to divide the 
inlet computational plane into 18 regions; 17 concentric annular regions 
forming the injection elements and holes and one ignition region. The 17 
injection regions correspond to the 17 rings of holes and elements and holes. 

Nine of these annular regions contain only hydrogen injection while in 
the remaining eight both hydrogen and oxygen injection occurs, and within each 
region fluid is injected nonuniformly. To most accurately model each region 
fluid should be injected with the proper momentum from computational nodes 
which form an area which is approximately the same as the effective area for 
that region. The remaining nodes in that region must have a zero velocity 
fixed boundary condition. The effective area is less than the actual combined 
areas of holes and elements. 
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Fig. 3-12 Fuel Prebumer Flow Rates Taken from Digital Transient Model (DTM) 



Fig. 3-13 Fuel Prebumer Flow Rates from DTM During Early Transient 
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Fig. 3-14 Fuel Prebumer Igniter Flow Rates from Digital Transient Model (DTM) 



Fig. 3-15 Fuel Preburner Igniter Flow Rates from DTM During Early Transient 
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Time from Start Command (sec) 


Fig. 3-16 Fuel Prebumer Average Chamber Pressure and Temperature 
During Start-Up Transient (Taken from DTM) 
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The mass flow rate, AM/At, through an area A of fluid having mass 
density of p and inlet speed V can be expressed as 


Am 

At 


p VA 


The ratio of dynamic pressure to pressure loss through the opening defines a 
discharge coefficient given by 


C d = 2 p v2/Ap 


Combining these two equations one expression for the product of A and 
is obtained, which corresponds to the effective area. 


AC d 


AM/At 

1/2 

0.67 (p Ap) 


In this equation the 0.67 includes unit conversion factors so that AC. is de- 

a 

termined in square inches. Data obtained from the SSME Power Balance Program, 
Revision F, provide a Ap for the entire face plate at FPL. Combining this 

with the appropriate mass flow rate at FPL one obtains an AC. for the entire 

2 d 
face plate of 3.32 in . This should be equal to the sum of the AC^ pro- 
ducts for all of the fuel holes and the annular injection elements. The actual 

2 

combined hole area is 0.792 in and that of the combined annular fuel elements 
2 

is 3.59 in . Assuming a C. of 1.0 for the holes we thus determine C. for 

d a 

the fuel elements to be 0.704. Because of the geometry of the oxygen injection 

path through the center of the injection elements it is reasonable to expect 

that the C. for these holes is essentially unity, 
d 


With this information an effective area which should be modeled by the 
computational grid for each annular inlet plane region can be compiled and is 
presented in Table 3-1. Note that if a mass flow rate is computed for each 
region and is input uniformly over that part of the grid the stream will have 
the incorrect momentum. It is absolutely necessary that the areas of the 
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Table 3-1 EFFECTIVE AREAS (IN SQUARE INCHES) OF EACH ANNULAR 
ZONE (360 DEGREES) OF PREBURNER INJECTOR FACE-PLATE 
FOR COMPUTATIONAL MODELING OF INLET PLANE 


Zone 

Fuel 

Elements 

Holes 

Baffles 

Oxidizer 
Injection Area 

1 


0.075745 



2 

0.54639 


0.0093 

0.335934 

3 


0.044639 



4 

0.48568 

0.041849 

0.0093 

0.298608 

5 


0.039060 



6 

0.42497 

0.036270 

0.0093 

0.261282 

7 


0.033479 



8 

0.36426 

0.030690 

0.0093 

0.223956 

9 


0.041780 



10 

0.30355 

0.037599 

0.0093 

0.186630 

11 


0.033420 



12 

0.24284 

0.029243 

0.0093 

0.149304 

13 


0.030750 



14 

0.18213 

0.025620 

0.00853 

0.111978 

15 


0.020500 



16 

0.12142 

0.015373 

0.00615 

0.074652 

17 


0.020500 

0.02214 
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regions of mass injection in each ring on the inlet plane closely match the 
appropriate effective area. Thus many of the inlet plane nodes in the grid 
will be "turned off." 

3.2.2 Chemistry Model 

During the startup transient it is not clear how much and where combustion 
or burning occurs. It cannot be assumed that during this time interval the 
chemistry is completed within a short distance of the face plate. Much 
research was performed to ascertain the most appropriate chemistry model to 
employ for computing a realistic flow-kinetic picture for this problem. 

It was determined that a non-equilibrium chemistry model be used to model 
the details of the kinetics. The eight reactive species are H 2 » C> 2 , H 2 0, 

H 2°2’ H ’ °* 0H ’ an< * H °2‘ A 20-reaction set was selected, eight of which apply 
to ignited hot portions of the flow field and 12 of which model ignition and 
cool flows. Tables 3-2 and 3-3 list these reactions along with their rate 
coefficients. Literature references are also included for completeness. 

A reduced set of reactions can be substituted for these 20-reactions for 
a less complete treatment. How much of the physics is lost, if any, can only 
be discovered by computing a test problem using both sets. This has not been 
done. 


The reduced set is obtained as follows: 

1. Use Reactions 1 through 6 in Table 3-2 

2. Use Reaction 9 in Table 3-3 

3. Rewrite Reaction 10 as 

H + 0 2 +M-»0.5H 2 +0 2 +M 

The eighth reaction is forward only, and the rate coefficient should be set 
equal to two thirds that for Reaction 10. Thus the minimum recommended 
chemistry treatment should contain eight reactions and six species. 
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Table 3-2 KINETIC REACTIONS (1 THROUGH 8) 


N Body Reaction 
(Reversible) 

Rate Coefficient 

(cm * particle - 1 i ^sec~ 1 units 

Ref. 

i.6h + h 2 -+h 2 o + h 

1.7 x 10“ 16 T'* 6 exp(- 3,300/RT) 

1,2 

2. II + O a - OH + O 

2.5 x 10- 7 r-°- 9 exi)(-17,330//?r) 

3,4 

3. 0 + II 2 -> OH + H 

1.8 x 10“ 2O T 2,8 exp(— 5, 920/ RT) 

4 

4. OH i OH — ► H 2 0 + O 

1.4 x 10- 2O r 2 - 7 exp(l,800/Rr) 

5, 6 

5. II + OH + M — ► H 2 0 + M 

2.4 x 10“ 26 T“ 2 

7,8 

G. II + II + M — H 2 + M 

2.o x io -30 t-' (M # H 2 ) 

9,2 


2.8 x 10“ 3I T“°- 6 (M = H 2 ) 

2 

7. H + O + M — OH + M 

4.0 x 10“ 29 T“ 1 

10 

8. O + O + M — 0 2 + M 

6.0 x 10~ 35 exp(l, 800/RT) 

11 


Table 3-3 KINETIC REACTIONS (9 THROUGH 20) 


N Body Reaction 
(Reversible) 

Rate Coefficient 

(cm 3 particle - l )( N ~ 1 ) sec - 1 units 

Ref. 

9. II a + 0 2 — OH + OH 

2.8 x 10“ 1 1 exp( — 48, 100/ RT) 

12 

10. H + 0 2 + M -> HOa + M 

2.0 x 10“ 3O T“° 8 

13 

11. II0 2 f OH — HaO + Oa 

1.6 x 10“ 11 cxp(l, 000/RT) 

14 

12. HOa + H-Ha + Oa 

4.2 x 10“ 1 1 exp(-700/RT) 

9,2 

13. II0 2 -b H — ► OH + OH 

1.9 x 10“ 10 exp(- 1, 000/RT) 

16 

14. HOa + O -> OH + 0 2 

3.3 x 10“ 11 cxp(400/ RT) 

17 

15. HOa + H0 2 — H 2 0 2 + 0 2 

3.2 x 10“ 12 cxp(2, 300/RT) 

15 

1G. HOa + H 2 — H 2 0 2 + H 

3.2 x 10“ 12 exp(-21, 000/RT) 

15 

17. HaOa + OH — H0 2 + H 2 () 

3.7 x 10“ 12 exp( — 500/RT) 

18 

18. HaOa + H — OH + H 2 0 

1.7 x 10“ u exp(-3, 600/RT) 

2 

19. HaOa + O — OH + H0 2 

8.7 x 10“ l3 exp(-3, 600/RT) 

20 

20. H 2 0 2 + M OH + OH + M 

1.2 x 10“ 7 cxp(-46, 300/RT) 

19 
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4. SUMMARY AND CONCLUSIONS 


A three-duct SSME Hot Gas Manifold geometry code has been developed for 
use by the computational mechanics staff of NASA-MSFC . This report describes 
the methodology of the program, makes recommendations on its implementation, 
and provides an input guide, and input deck listing, and a source code 
listing. The code listing is strewn with an abundance of comments to assist 
the user in following its development and logic. A working source deck will 
be provided to NASA-MSFC on the EADS network upon request. 

A thorough analysis has been made of the proper boundary conditions and 
chemistry kinetics necessary for an accurate computational analysis of the 
flow environment in the SSME fuel side prebumer chamber during the initial 
startup transient. Pertinent results have been presented to facilitate in- 
corporation of these findings into an appropriate CFD code. The computation 
must be a turbulent computation, since the flow field turbulent mixing will 
have a profound effect on the chemistry. Because of the additional equations 
demanded by the chemistry model it is recommended that for expediency a simple 
algebraic mixing length model be adopted. 

Performing this computation for all or selected time intervals of the 
startup time will require an abundance of computer CPU time regardless of the 
specific CFD code selected. 
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Appendix A 


INTRODUCTION 


The geometry input guide is presented in two sections: (1) a definition 
of terminology commonly used for inputting and describing the geometry, and (2) 
a summary of card types used to input the geometry and a detailed description 
of the associated parameters and their input values. 

We begin with an overview of how to apply the program. The flowfield 
domain is divided into zones in order to simplify the input necessary to 
describe the complicated geometry. Each zone contains its own coordinate 
system. A zone consists of at least one section. A multiple section zone is 
used to describe Zones 4 and 5. Sections may be added in only one internal 
coordinate direction which is called the marching direction since this is the 
direction in which the grid is created. Each section is described using 
points, edges, and surfaces. An edge may consist of from one to ten segments. 

A segment or a surface may require additional input depending on its type. 

The second section presents a detailed description of the input param- 
eters. Each card type is listed in the order of input with its associated 
parameters. Each parameter is identified as to its usage in the program with 
the options of each shown. Reference to Fig. A-l or Table A-l may be necessary 
to explain some of the input parameters and their order of input. All of the 
card types are not necessarily input for a specific zone. 

Card type 9 may be used when other than an equal distribution of nodes is 
desired, etc. Whereas card type 13 is necessary if there will be more than one 
segment per edge. And cards type 10 and 11 are used if additional information 
is needed to describe a segment or a surface. Certain of the input parameters 
on early cards dictate which of the later cards are read in. 
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DEFINITION 

Edge 

Map 

Node 

Point 

Section 

Segment 

Surface 

Zone 

D21 13 


OF TERMINOLOGY 


An edge consists of from one to five segments. Four edges are 
used to describe a 2D section. Twelve edges are used to describe 
a 3D section. 

The geometry maps a point from n space into real space. When 
describing a surface mapping we could say map = 2 refers to a 
planar n space surface being mapped onto a cylindrical surface 
in real space. 

At each intersection of t) coordinates a node is generated by 
the program forming the grid which will describe the flowfield 
domain. 

The comers of a section are called "points." The location and 
initial flow directions are input for each point. There are 
eight points. 

A zone may be further subdivided into sections. Each zone 
consists of at least one section. A zone may only be subdivided 
into sections in the direction in which the grid is generated 
(marching direction). Points, edges or a surface in common with 
a previous section within a zone are not reinput. Sectioning is 
done primarily to reduce the amount of input. Each section is 
described using points, edges, and surfaces. 

An edge is subdivided into as many as five segments. A segment 
may be a straight line, a circular arc, a helical coil, a 
trigonometic function of angle or length, a cubic spline, or user 
defined. 

A three-dimensional section will consist of six surfaces which 
form a generalized hexahedron. A surface may be planar, 
cylindrical, an edge of revolution, or user defined. 

The flowfield domain may be subdivided into zones. Zones are 
generated independently and are the fundamental building block. 
Each zone contains its own n coordinate system. 

Localized coordinate directions within a zone. These co- 
ordinates describe a cube in n space. 0 < nj < 1 
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SUMMARY OF CARDS 


Card Type 


Parameter List/Format 


Flowf ield 
Parameters 

1 ITITLE(I) , 1=1,80 

(20A4) 

2 NZONE , MAPTEN, INCHES 

(815) 

3 IWRTI, IWRTC , IWRTN 

(315) 


Zone Parameter 


4 NSECT 

(15) 

Section 

Parameters 

5 MAPEDGE(I) , 1=1,12 

(1215) or (3(4110)) 

6 MAPSIDE(I) , 1=1,6 

(615) 

7 MARCH 

(15) 

8 (NMBRNDS(I) , 1=1,3), (ISTRTCH(I), 1=1,3) 

(615) 

9 STRETCH(I) , 1=1,3 

(3E10.4) 

10 [ ( COEFE ( I , K , J ) , 1=1,8), K=l,5] , J=1,4(2D) OR 12(3D) 

(8E10.4) 

11 [COEFS(I, J) ,1=1,8] , J=l,6 

(8E10.4) 

12 (POINT ( I, J) , 8=15), J=l, 8 

(8E10.4) 

13 [ ( SEGMAZ ( I , K , J ) , 1=1,5), ETAMAX(K.J), K=l,4] J=l, 12 

(6E10.4) 
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CARD AND VARIABLE DESCRIPTIONS 


Flowfield Parameters 


CARD TYPE 1 Problem Identification Label Format(20A4) 

ITITLE Alphanumeric information used for identifying the flowfield 
geometry. Columns 1-80 are read and printed only. 


CARD TYPE 2 Problem Option Controls Flags Format(5I5) 

NZONE The number of zones into which the flowfield geometry is 
divided. The maximum number of zones is 99. 

HAPTEN This option determines the maximum number of segments which 
will be input per edge. 

= 0 Five segments per edge Format(12I5) 

= 1 Ten segments per edge Format(3(4I10) ) 

INCHES This option specifies the dimensions of the coordinates being 
input. The output data will be written in feet for compati- 
bility with the INTEGRATION program. 

= 0 Dimensions in feet 
= 1 Dimensions in inches 


CARD TYPE 3 Output Print Options Format(3I5) 

IWRTI This option may be used when either grid plots or nodal printout 
are insufficient to determine the cause of an error in the 
geometry. Intermediate debug printout may be generated in 
subroutines EDGE and SURFACE. Subroutine EDGE prints the 
following information: node number, edge number, position on 
the edge, n value at this position on the edge, a ratio which 
gives the relative location on a segment, and the maximum n 
value for this segment of the edge Subroutine SURFACE prints 
the following information: node number, surface number, and 
position on the surface. It also prints the following informa- 
tion for an edge of revolution: node number, surface number, 
position on the edge of revolution to be revolved, a vector 
tangent to the edge at this position, edge number, the distance 
along the axis of revolution from the origin of the relative 
coordinate system to a perpendicular from the axis to the point 
on the revolved edge, and the radius along this perpendicular 
from the axis to the revolved point on the edge. IWRTI should 
be chosen carefully to restrict the amount of printout. 
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=0 No debug print 

= N > 0 Intermediate error printout for every N fc ^ node 
= N > Total number of nodes - Prints the number of nodes 
stored plane. 

IWRTC This option controls the printout of connectivities for all 

nodes specified for MATING. Incorrect connectivities result 
in various types of errors in the INTEGRATION program and a 
person familiar with it should be consulted if this type of 
error is suspected. 

=0 Do not print mated connectivities 
= 1 Print mated connectivities 

IWRTN This option controls the amount of nodal printout for a run. 

The nodal information will be printed will be node number and 
position. For large problems with many thousands of nodes, 
IWRTN should be chosen carefully to restrict the amount of 
paper produced. A plot of the grid is usually more instructive 
than a massive nodal printout for locating errors in the code. 

=0 No printout 

= N > 0 Print nodal point information every N th node 


Zone Parameter 


CARD TYPE 4 Number of Sections Format (15) 

NSECT The number of sections in this zone. A zone may only be 
sectioned in the n direction which corresponds to MARCH. 
There is no limit to the number of sections within a zone. 


c 

Section Parameters 


CARD TYPE 5 Edge Shape Function Indicators Format (1215) or 3(110) 

MAPEDGE ( I ) . 1=1. 12 

These are packed integer flags that specify which edge shape 
functions will be used for the current section. The edges are 
input in numerical order. The edge numbers are defined 
according to Fig. A-l. The user should study this figure 
before inputting the geometry. 
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Each of the edges may consist of up to ten segments with each 
of these segments having its own shape function. The value of 
MAPEDGE(I) can consist of up to ten integers packed into one 
word MAPEDGE(I). MAPTEN specifies the maximum number of seg- 
ments per edge. The edge shape function indicators for each 
segment are input in chronological order of increasing n for 
each edge with the final packed integer being right adjusted. 

For example, if MAPEDGE(4) = 112, then edge 4 consists of three 
segments: the first segment is type 1; the second segment is 
type 1; and the third segment is type 2. If only one segment 
describes an edge, then only one indicator is used, right 
adjusted. 

A library of edge shape functions indicators for the HGM 
GEOMETRY program follow. If any edge shape function other than 
a linear segment is specified, then edge coefficients (COEFE(I)) 
must be input. Card type 1 is used to define the analytical 
function describing a segment. 

= 1 Linear segment 
= 2 Circular arc 
= 3 Edge of revolution 
= 4 Special segment 
= 5 Special segment 


CARD TYPE 6 Surface Shape Function Indicators Format(6I5) 

MAPSIDE(I) .1=1,6 

These are integer flags that specify which surface shape 
functions will be used for the current section. These flags 
are input only for three-dimensional problems since two- 
dimensional geometries are defined completely by the edge 
functions. The surfaces are input in numerical order. The 
surface numbers are defined in Fig. A-l. The user should study 
this figure before inputting the geometry. An edge of revolu- 
tion requires the input of surface coefficients (COEFS(I)) on 
card type 15 to define a relative origin on the axis, the axis 
of revolution, and the direction of revolution. 

= 1 Planar surface 
= 2 Cylindrical surface 

= 3 Special surface (user defined) 

= 4 Edge of revolution (input COEFS(I)) 

= 5 Holes in bowl surface 
v =6 Hole side of duct surface 
= 7 Duct outer surface 


(input COEFE(I) ) 
(input COEFE(I) ) 
(input COEFE(I) ) 
(input COEFE(I) ) 
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CARD TYPE 7 Node Numbering Sequence Specs Format(6I5) 

MARCH (15) The value of MARCH determines the node generation and hence the 
node numbering sequence, (default =1) 

The numbering sequence corresponding to follows. 

= 1 TI3 , t» 2 . Hi (default) 

= 2 m, TI3, r\ 2 

= 3 Ti2 » 'll* ’ll 


CARD TYPE 8 Node Distribution Parameters 


Format(6I5) 


NMBRNDS(I) .1=1 3 

NMBRNDS(I) is the number of nodes in the Tij direction for 
the current section. The limit is 200 nodes in any coordinate 
direction. This may be changed in the program by respecifying 
the ETAS(3,200) array. 


ISTRTCH(I) , 1=1 3 

This option gives the user control over the node distribution 
in each of the coordinate directions. 


= 0 
= 1 

= 2 

= 3 

= 4 

= 5 

= 6 

= 7 


Uniform spacing 

Input actual nj values for NMBRNDS(I) nodes 
(input ETAS ( I ) ) 

Decrease spacing in nj direction. Input a stretching 
factor greater than 0.0 in STRETCH(I). 

Increase spacing in rij direction. Input a stretching 
factor greater than 0.0 in STRETCH(I). 

Double stretching. Input a stretching factor greater than 
0.0 in STRETCH(I) . Use an odd number of nodes. 

Decrease spacing in Hi direction. Input minimum grid 
spacing as a percentage of the total length in STRETCH(I). 
Increase spacing in direction. Input minimum grid 
spacing as a percentage of the total length in STRETCH(I). 
Double stretching. Input minimum grid spacing as a 
percentage of the total length in STRETCH(I). Use an odd 
number of nodes. 


If ISTRTCH(I) = 1, input a set of cards type 13 for each n 

direction to be input. 

If ISTRTCH(I) >2, input card type 12. 
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CARD TYPE 9 Option for Stretching Function Format (3E10. 4) 

(input when ISTRTCH(I) > 2) 

STRETCH(I) .1=1. 3 

This parameter is input for each coordinate direction 
designated for stretching by ISTRTCH(I) > 2. 

Example: Several stretching functions will be demonstrated using 21 

points for comparison. Note, that total length = 10.0 for 
ISTRTCH = 6 and 7 . 

ISTRTCH = 3 STRETCH =2.0 



STRETCH =4.0 



STRETCH =6.0 


STRETCH =8.0 


STRETCH =10.0 


STRETCH =2.0 


1 


I 


[ 


ISTRTCH = 4 
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STRETCH =4.0 


STRETCH = 


STRETCH =8.0 


ISTRTCH « 6 

ML 




STRETCH =10.0 


STRETCH 


.002 


STRETCH 


.004 


A-10 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 






LMSC-HEC TR F225755 


STRETCH = .010 


ISTRTCH 


7 


STRETCH = .002 




STRETCH a .008 


STRETCH = .010 



CARD TYPE 14 Coefficients for Edge Shape Functions Format (8E10. 4) 

(input for each segment in MAPEDGE(I) > 1) 

COEFE(I) .1=1.8 

These coefficients are used to describe the edge shape 
functions for each segment of the current section. The 
coefficients for each segment are input on separate cards in 
the same order as the indicators on card type 5. 
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Input Parameters 

COEFE(I) , 1=1,3 are the x, y, and z coordinates of the 
center of the arc. 

COEFE(I) , 1=1,3 are coordinates of the center 
of the arc. 

COEFE(I) , 1=4,6 are components of the unit vector 
along the axis. Direction according to right hand 
rule when a vector is revolved from point 1 to point 2. 


CARD TYPE 15 Coefficients for Surface Shape Functions Format (8E10. 4) 
(input when MAPSIDE(I) = 4 and IDIM = 3) 

COEFS(I) .1=1.8 

These are the coefficients defining the surface shape functions 
for each surface formed by an edge of revolution in the current 
section. Each surface which has MAPSIDE(I) = 4 on card type 8 
is input on a separate card in the same order as they occur on 


Input Parameters 

Surface formed by revolving an edge about an 
axis . 

COEFS(I), 1=1,3 a point on the axis of revolu- 
tion which becomes the origin of a local co- 
ordinate system. This point must lie outside 
of the projection of the edge onto the axis of 
revo lution. 

COEFS(I), 1=4,6 are components of the unit 
vector along the axis of revolution in the 
direction of increasing n. 

COEFS( 7 ) indicates the tij direction in 
which the edge will revolve. 


CARD TYPE 16 Coordinates of Points and flow direction Format (5E10. 4) 

POINT(I.J). J = Point number 

These parameters are the coordinates and flow direction at each 
corner of a general hexahedral(3D) . Figure A-l shows this con- 
figuration with the points numbered from 1 to 8. There are 
eight cards of type 16 to be input. 


card type 8. 

Map Type 

4 Surface of 
Revolution 


Map Type 

2 Circular Arc 

3 Edge of 
Revolution 
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P0INT(1,J) - the x coordinate of point J 

POINT(2,J) - the y coordinate of point J 

P0INT(3,J) - the z coordinate of point J 

POINT(4 , J) - the flow angle 0 at point J 

POINT(5,J) - the flow angle $ at point J 

Important Note: All cards type 12 are not input consecutively. 

They are grouped with cards type 13. See Table A-l for the 
exact sequence of card types 11 and 13. 


CARD TYPE 17 Segment Extremals for Edges 


Format(6E10.4) 


SEGMAX(I ,K. J) . K = Segment Number. J = Edge Number 

Each edge may be segmented up to five times. Therefore, cards 
type 13 are repeated for each successive segment of edge J. 

Each segment must be input on a separate card type 13. The 
extremal for the final segment of an edge is not to be input 
since this point is already defined by the POINT(I) input. The 
number of cards type 13 for each edge will thus be one less 
than the number of segments on that edge. In particular, if an 
edge consists of only one segment, no cards of type 13 are 
input for that edge. 


See Table 1 for the input order of card types 12 and 13. Each 
POINT(I) is input on a single card, followed by up to five 
cards containing the extremals. 


SEGMAX(1,K, J) 
SEGMAX(2,K,J) 
SEGMAX(3,K,J) 
SEGMAX(4 ,K, J) 
SEGMAX(5 , K, J) 


The extremal x coordinate 
segment of edge J. 

The extremal y coordinate 
segment of edge J. 

The extremal z coordinate 
segment of edge J. 

The extremal flow angle 0 
segment of edge J. 

The extremal flow angle $ 
segment of edge J. 


for 

the 

K th 

for 

the 

K th 

for 

the 

K th 

for 

the 

K th 

for 

the 

K th 


ETAHAX(K.J), K= Segment Number, J= Edge Number 

The maximum value of the nj coordinate on the segment 
of edge J. Input a negative value when defining a fold line. 
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Table A-l INPUT SEQUENCE FOR POINTS AND EXTREMALS OF SEGMENTS: 
CARDS TYPE 16 AND 17 


Order 

Card Type 

Description 






1 

12 

Point 1 







2 

13 

Extremals 

for 

each 

segment 

of 

Edge 

1 

3 

12 

Point 2 







4 

13 

Extremals 

for 

each 

segment 

of 

Edge 

2 

5 

12 

Point 3 







6 

13 

Extremals 

for 

each 

segment 

of 

Edge 

3 

7 

12 

Point 4 







8 

13 

Extremals 

for 

each 

segment 

of 

Edge 

4 

9 

12 

Point 5 







10 

13 

Extremals 

for 

each 

segment 

of 

Edge 

5 

11 

12 

Point 6 







12 

13 

Extremals 

for 

each 

segment 

of 

Edge 

6 

13 

12 

Point 7 







14 

13 

Extremals 

for 

each 

segment 

of 

Edge 

7 

15 

12 

Point 8 







16 

13 

Extremals 

for 

each 

segment 

of 

Edge 

8 

17 

13 

Extremals 

for 

each 

segment 

of 

Edge 

9 

18 

13 

Extremals 

for 

each 

segment 

of 

Edge 

10 

19 

13 

Extremals 

for 

each 

segment 

of 

Edge 

11 

20 

13 

Extremals 

for 

each 

segment 

of 

Edge 

12 
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Notes on Input Sequence of Cards Type 4 Through 13 


• Card types 5-13 are repeated for each section in a zone. 
Sections can be added only in the direction of MARCH. 

• Card type 4 is repeated for each zone. 

• Cards type 9-13 are not re-input for points, edges and 
surfaces in common with the preceding section. Since these 
have already been input, the code transfers the entries 
from one section to another. 

• The first section of each new zone, however, must be input 
since zones are considered to be independent. 
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Appendix B 

HGM GRID CODE INPUT LISTING 
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Appendix B 


ZONE 1 ( TAD ) FOR 3 -DUCT HGM GEOMETRY 


1 

3 

1 

1 






10000 

0 

0 







1 

1311 


1 1311 


1 



ZONE 


3 


3 3 


3 





1311 


1 1311 


1 




1 

4 

1 

4 4 4 






3 









94 

21 

72 

0 4 0 






0.0 


8.0 

0.0 






1.35 


0.0 

- 6.0 

0.0 


- 1.0 

0.0 


1.35 


0.0 

- 5.875 

0.0 


- 1.0 

0.0 


3.31 


0.0 

0.0 

1.0 


0.0 

0.0 


6.07 


0.0 

0.0 

1.0 


0.0. 

0.0 


6.07 


0.0 

0.0 

1.0 


0.0 

0.0 


3.31 


0.0 

0.0 

1.0 


0.0 

0.0 


1.35 


0.0 

6.0 

0.0 


1.0 

0.0 


1.35 


0.0 

5.875 

0.0 


1.0 

0.0 


6.07 


0.0 

0.0 

1.0 


0.0 

0.0 

3.0 

6.07 


0.0 

0.0 

1.0 


0.0 

0.0 

3.0 

6.07 


0.0 

0.0 

1.0 


0.0 

0.0 

3.0 

6.07 


0.0 

0.0 

1.0 


0.0 

0.0 

3.0 

3.31 


0.0 

- 5.55 






1.35 


0.0 

- 5.55 

18.0 





1.35 


0.0 

- 6.45 

56.0 





4.02 


0.0 

- 6.5 

81.0 





6.07 


0.0 

- 6.5 






6.07 


0.0 

- 7.434 






1.35 


0.0 

- 4.525 

18.0 





1.35 


0.0 

- 7.225 

55.0 





4.02 


0.0 

- 7.434 

81.0 





3.31 


0.0 

- 4.525 






3.31 


0.0 

5.55 






6.07 


0.0 

6.5 






6.07 


0.0 

7.434 






3.31 


0.0 

4.525 






1.35 


0.0 

5.55 

18.0 





1.35 


0.0 

6.45 

56.0 





4.02 


0.0 

6.5 

81.0 





1.35 


0.0 a 

4.525 

18.0 





1.35 


0.0 

7.225 

55.0 





4.02 


0.0 

7.434 

81.0 






ZONE 2 ( 1/2 

BOWL ) 

FOR 

3 -DUCT HGM 

GEOMETRY 

1 3 

1 

1 



10000 0 

0 




1 




ZONE 

111111111 


1 

121111 

1 

3 


3 

3 

3 

111111111 


1 

121111 

1 

1 4 

1 

4 

4 4 



3 
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59 21 

36 

0 4 0 

0.0 

8.0 

0.0 

9.445 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

9.445 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

. 0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

6.07 

0.0 

- 6.5 

6.8015659 

0.0 

- 6.5 

7.8 

0.0 

- 6.5 

9.1 

0.0 

- 6.62 

9.78 

0.0 

- 6.62 

13.555 

0.0 

- 6.37 

13.8395195 

0.0 

- 6.2211618 

13.8395195 

0.0 

- 6.5734119 

6.8015659 

0.0 

- 7.434 

13.39 

0.0 

- 6.8329404 

6.07 

0.0 

- 7.434 

6.07 

6.5 

0.0 

13.8395195 

6.2211618 0.0 

13.8395195 

6.5734119 0.0 

6.07 

7.434 

0.0 

6.8015659 

6.5 

0.0 

7.8 

6.5 

0.0 

9.1 

6.62 

0.0 

9.78 

6.62 

0.0 

13.555 

6.37 

0.0 

6.8015659 

7.434 

0.0 

13.39 

6.8329404 0.0 


1.0 

0.0 

0.0 


1.0 

0.0 

0.0 


1.0 

0.0 

0.0 


1.0 

0.0 

0.0 


1.0 

0.0 

0.0 

3.0 

1.0 

0.0 

0.0 

3.0 

1.0 

0.0 

0.0 

3.0 

1.0 

0.0 

0.0 

3.0 

4.0 




11.0 




21.0 




26.0 




54.0 





4.0 

53.0 


4.0 

11.0 

21.0 

26.0 

54.0 
4.0 

53.0 


ZONE 3 ( 1/2 BOWL ) FOR 3 -DUCT HGM GEOMETRY 
13 11 

10000 0 0 
1 


111111111 


1 

121111 


1 



3 


3 


3 


3 



111111 


1 


121 


1 



1 4 

1 

5 

4 

4 





3 0 

1 








4.0 


20.0 


46.0 





13.0 


41.0 


58.0 





3.0 


4.0 


4.0 





3.0 


4.0 


4.0 


3.0 



59 21 

72 

0 

4 

0 





0.0 

8.0 


0.0 






9.445 

0.0 


0.0 






0.0 

0.0 


0.0 


1.0 


0.0 

0.0 

0.0 

0.0 


0.0 


1.0 


0.0 

0.0 

0.0 

0.0 


0.0 


1.0 


0.0 

0.0 

0.0 

0.0 


0.0 


1.0 


0.0 

0.0 

9.445 

0.0 


0.0 







ZONE 3 
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0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

3.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

3.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

3.0 

0.0 

0.0 

0.0 

1.0 

0.0 

0.0 

3.0 

6.07 

6.5 

0.0 





6.8015659 

6.5 

0.0 

3.0 




7.8 

6.5 

0.0 

11.0 




9.1 

6.62 

0.0 

21.0 




9.78 

6.62 

0.0 

26.0 




13.555 

6.37 

0.0 

55.0 




13.839518 

6.221162 

0.0 





13.839518 

6.573411 

0.0 





6.8015659 

7.434 

0.0 

4.0 




13.39 

6.8329404 

0.0 

53.0 




6.07 

7.434 

0.0 





6.07 

0.0 

6.5 





13.839518 

0.0 

6.221162 





13.839518 

0.0 

6.573411 





6.07 

0.0 

7.434 





6.8015659 

0.0 

6.5 

4.0 




7.8 

0.0 

6.5 

11.0 




9.1 

0.0 

6.62 

21.0 




9.78 

0.0 

6.62 

26.0 




13.555 

0.0 

6.37 

55.0 




6.8015659 

0.0 

7.434 

4.0 




13.39 

0.0 

6.8329404 

53.0 





ZONE 4 (TOP T - D ) FOR 3 -DUCT HGM GEOMETRY 


1 

3 1 

1 






10000 

0 0 







3 







ZONE 4 


4 

5 


3 

5 




4 

4 


3 

3 




4 

5 


3 

5 



7 

6 7 

1 

7 

7 




2 








27 20 29 

4 

0 

4 




7.0 

0.0 


7.0 





0.0 

0.0 


0.0 

0.0 

0.0 

0.0 


0.0 

0.0 


0.0 

0.0 

0.0 

0.0 


9.445 

5.564 


0.0 

- 0.25506 

- 0.092215 

- 0.96252 


0.0 

0.0 


0.0 

0.0 

0.0 

0.0 


0.0 

0.0 


0.0 

0.0 

0.0 

0.0 


0.0 

0.0 


0.0 

0.0 

0.0 

0.0 


9.445 

5.564 


0.0 

- 0.25506 

- 0.092215 

- 0.96252 


9.445 

5.564 


0.0 

0.25506 

0.092215 

0.96252 


0.0 

0.0 


0.0 

0.0 

0.0 

0.0 


0.0 

0.0 


0.0 

0.0 

0.0 

0.0 


9.445 

5.564 


0.0 

0.25506 

0.092215 

0.96252 


0.0 

0.0 


0.0 

0.0 

0.0 

0.0 


. 0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

8.691001 7.428427 

2.549918 
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11.362797 


7.428428 

1.841919 







12.66737 


7.79672 

7.24362 







10.23170 


7.79672 

7.88905 







9.568123 


4.613040 

6.399753 







12.403425 


4.49657 

5.768945 







12.73830 


4.80664 

7.51127 







10.17804 


4.92167 

8.17869 








3 

1 

3 


1 





3 

3 

3 


3 





3 

1 

3 


1 




4 

-i 

l 

4 

1 4 

4 






27 13 

29 

4 0 

4 






7.0 


0.0 

7.0 







9.445 


5.564 

0.0 


-0.25506 


-0.092215 

-0.96252 


9.445 


5.564 

0.0 


-0.25506 


-0.092215 

-0.96252 


9.445 


5.564 

0.0 


0.25506 


0.092215 

0.96252 


9.445 


5.564 

0.0 


0.25506 


0.092215 

0.96252 


9.445 


5.564 

0.0 


-0.25506 


-0.092215 

-0.96252 

1.0 

9.445 


5.564 

0.0 


0.25506 


0.092215 

0.96252 

1.0 

9.445 


5.564 

0.0 


0.25506 


0.092215 

0.96252 

3.0 

9.445 


5.564 

0.0 


-0.25506 


-0.092215 

-0.96252 

3.0 

14.18025 


8.61610 

11.69913 







11.12323 


8.61610 

12.50921 







14.26927 


4.86326 

12.03505 







11.05588 


5.00762 

12.87274 








3 

1 

3 


1 





3 

3 

3 


3 





3 

1 

3 


1 




4 

1 

4 

1 4 

4 






27 13 

29 

4 0 

4 






7.0 


0.0 

7.0 







9.445 _ 


5.564 

0.0 


-0.25506 


-0.092215 

-0.96252 


9.445 


5.564 

0.0 


-0.25506 


-0.092215 

-0.96252 


9.445 


5.564 

0.0 


0.25506 


0.092215 

0.96252 


9.445 


5.564 

0.0 


0.25506 


0.092215 

0.96252 


9.445 


5.564 

0.0 


-0.25506 


-0.092215 

-0.96252 

1.0 

9.445 


5.564 

0.0 


0.25506 


0.092215 

0.96252 

1.0 

9.445 


5.564 

0.0 


0.25506 


0.092215 

0.96252 

3.0 

9.445 


5.564 

0.0 


-0.25506 


-0.092215 

-0.96252 

3.0 

15.51900 


9.10020 

16.75236 







12.46198 


9.10020 

17.56244 







15.60802 


5.34736 

17.08828 







12.39463 


5.49172 

17.92597 








ZONE 5 
1 
0 
3 


7 

2 

27 


(MID T-D) FOR 3-DUCT HGM GEOMETRY 
3 1 1 

0 0 


4 5 3 
4 4 3 
4 5 I 
6 117 7 


5 

3 

5 


19 16 4 0 3 
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7.0 

0.0 

5.0 





0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


7.32766 

0.0 

0.37334 

- 0.40139 

0.0 

- 0.91591 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


7.32766 

0.0 

0.37334 

- 0.40139 

0.0 

- 0.91591 


7.32766 

0.0 

0.37334 

0.40139 

0.0 

0.91591 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

9.040005 

1.695449 

7.695033 





11.889400 

1.695342 

7.307725 





12.3459 

1.5014 

8.6854 





12.5633 

1.5014 

9.6966 





8.233353 

0.0 

7.796413 





12.667614 

0.0 

7.201867 





12.98742 

0.0 

8.40500 





9.39706 

0.0 

9.97844 






3 13 1 


3 3 3 3 

1111 
4 1114 4 


2 


27 25 

16 

4 0 

3 




7.0 

0.0 

5.0 





7.32766 

0.0 

0.37334 

- 0.40139 

0.0 

- 0.91591 


7.32766 

0.0 

0.37334 

0.40139 

0.0 

0.91591 


7.32766 

0.0 

0.37334 

0.40139 

0.0 

0.91591 


7.32766 

0.0 

0.37334 

- 0.40139 

0.0 

- 0.91591 

1.0 

7.32766 

0.0 

0.37334 

0.40139 

. 0.0 

0.91591 

3.0 

7.32766 

0.0 

0.37334 

- 0.40139 

0.0 

- 0.91591 

3.0 

14.7473 

1.8845 

13.3654 





11.8507 

1.8845 

14.6346 





15.55268 

0.0 

13.01283 





11.04640 

0.0 

14.98767 






3 13 1 


3 3 3 3 

1111 
4 1114 4 


2 


27 

4 16 

4 0 

3 




7.0 

0.0 

5.0 





7.32766 

0.0 

0.37334 

- 0.40139 

0.0 

- 0.91591 


7.32766 

0.0 

0.37334 

0.40139 

0.0 

0.91591 


7.32766 

0.0 

0.37334 

0.40139 

0.0 

0.91591 


7.32766 

0.0 

0.37334 

- 0.40139 

0.0 

- 0.91591 

1.0 

7.32766 

0.0 

0.37334 

0.40139 

0.0 

0.91591 

3.0 

7.32766 

0.0 

0.37334 

- 0.40139 

0.0 

- 0.91591 

3.0 

14.9963 

1.8845 

13.9344 





12.0997 

1.8845 

15.2036 





15.80235 

0.0 

13.58253 





11.29607 

0.0 

15.55737 
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Appendix C 


y 

PROGRAM HGMGEOM( INPUT , OUTPUT, FILE19 , FI LE20 , TAPES ,TAPE6=0UTPUT , 
& TAPE 19 = FI LEI 9 , TAPE 20= FILE 20 ) 

PROGRAM HGMGEOM 


DEVELOPED BY THE COMPUTATIONAL MECHANICS SECTION 
LOCKHEED ENGINEERING CENTER, HUNTSVILLE, ALABAMA. 

L.A. NICHOLSON 


TAPE19 GEOMETRY SCRATCH FILE 
TAPE20 GEOMETRY DATA 


COMMON /COUNTER/ 
COMMON /INITA/ 
COMMON /ZONING/ 
COMMON /UNITS/ 


NODESAV, NODETOT , NBNODES , NPLANE 

IDIM,MAPTEN, INCHES 

ISECT , NSECT , IZONE , NMBRNDS ( 3 ) 

NU5,NU6,NU19,NU20,NU21 


INITIALIZE PROGRAM 


CALL INITIAL(NZONE) 


ZONE 


DO 300 IZONE=l , NZONE 


READ( NU5 , 1000 ) NSECT 


SECTION 

DO 200 ISECT=1, NSECT 

READ INPUT FOR EACH SECTION 

CALL INPUT 

GENERATE GRID FOR EACH SECTION 

200 CALL GRID 

CALL REWRITE 
300 CONTINUE 

STOP 

—FORMAT STATEMENTS 
1000 FORMAT ( 15) 

END 
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SUBROUTINE INITIAL(NZONE) 


COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


/COEFF/ COEFE ( 8,10,12) ,COEFS ( 8 , 6 ) , NMBRSEG( 1 2 ) 
/COUNTER/ NODESAV, NODETOT , NBNODES , NPLANE 
/HEADER/ ITITLE( 20) ,LINE 
/INITA/ IDIM,MAPTEN, INCHES 

/INITB/ IBOX (4,6) , IRAY(4,3) , NRAY , IPT1 ( 12 ) ,IPT2(12) 
/INITC/ PI , RADDEG 
/IWRITE/ IWRTC, IWRTI, IWRTN 
/NODNMBR/ NODENUM( 10000) # TOL(3) 

/UNITS/ NU5 , NU6 , NU19 , NU20 , NU21 


DIMENSION NAME( 20 ) 

DATA PI /3. 141592654/ 

DATA RADDEG /57 . 29577951/ 

DATA IPT1 /l, 2, 4, 1,1, 2, 3, 4,5, 6, 8, 5/ 

DATA IPT2 /2, 3, 3, 4, 5, 6, 7, 8, 6, 7, 7, 8/ 

DATA I RAT / 1 , 3 , 9 , 1 1 , 4 , 2 , 1 2 , 10 , 5 , 6 , 8 , 7/ 

DATA IBOX /4, 1,2, 3, 5, 1,6, 9, 12, 9,10, 11, 8, 3, 7, 11, 5, 4, 8, 12, 6, 2, 


1 2 3 4 5 6 7 8 9 10 11 12 


IPT1 =124112345685 
IPT2 = 233456786778 

IRAY = 14 5 

3 2 6 

9 12 8 

11 10 7 


IBOX 


4 5 12 8 5 6 

1 1 9 3 4 2 

2 6 10 7 8 7 

3 9 11 11 12 10 


— DYNAMIC DIMENSIONER INPUT 


— INITIALIZE 

NU5 = 5 
NU6 = 6 
NU19 = 19 
NU20 = 20 
NU21 = 21 
NPLANE = 0 
NBNODES = 0 
NINODES = 0 
NODETOT = 0 
NODESAV = 0 
C 

DO 40 1=1,8 
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DO 30 J=1 , 6 
30 C0EFS(I,J) = 0.0 

DO 40 J=l,10 
DO 40 K=l,12 
40 COEFE ( I , J , K) = 0.0 


PRINT HEADER 

WRITE(NU6, 1120) 


NZONE 

IDIM 

MATING 

IWRTI 

IWRTC 

IWRTN 

MAPTEN 

INCHES 


NUMBER OF ZONES TO BE USED TO GENERATE GRID 

2 FOR 2D 

3 FOR 3D 

0 NO ZONE MATING 
0 NO INTERMEDIATE PRINT 
0 NO CONNECTIVITY PRINTOUT 
0 PRINT ONLY BOUNDARY NODES 
N PRINT EVERY NTH INTERIOR NODE 

0 FIVE SEGMENTS PER EDGE 

1 TEN SEGMENTS PER EDGE 

0 COORDINATES INPUT IN FEET 

1 COORDINATES INPUT IN INCHES 


READ( NU5 , 1000 ) ITITLE 
READ(NU5, 1010) NZONE, IDIM, 

& MAPTEN, INCHES 


READ( NU5 ,1010) IWRT I , IWRTC , IWRTN 


NRAY = 2* ( IDIM - 1) 

WRITE (NU6, 11 30) ITITLE 
WRITE(NU6, 1140) NZONE, IDIM, 

S MAPTEN, INCHES 

REWIND( NU20 ) 

DRAW PICTURE 

CALL PICTURE ( IDIM) 

RETURN 

FORMAT STATEMENTS 

1000 FORMAT ( 20 A4 ) 

1010 FORMAT ( 1615) 

1120 FORMAT ( 1H1 

3 / 40X, 34H LOCKHEED-HUNTSVILLE 

4 / 40X,34H VAX 11/785 VERSION 

5 / 40X, 32H HGM GEOMETRY ) 

1130 FORMAT { // 5X,20A4 ) 

1140 FORMAT (// 21H GEOMETRY PARAMETERS: 

1 // 45H NZONE IDIM 

2 MAPTEN INCHES /, 8110 ) 
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C 

END 

SUBROUTINE INPUT 


COMMON 

COMMON 

COMMON 

COMMON 

& 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

& 

COMMON 

COMMON 


/COEFF/ COEFE(8,10, 12) ,COEFS(8,6) , NMBRSEGf 12) 

/COUNTER/ NODESAV, NODETOT, NBNODES , NPLANE 
/HEADER/ ITITLE< 20) ,LINE 

/HOTGAS/ IHGM , ETA1 1 , ETA12 , ETA1 3 , ETA31 , ETA3 2 , ETA33 , 
STR1 1 , STR12 , STRl 3 , STR31 , STR3 2 , STR33 , STR34 
/INITA/ IDIM,MAPTEN, INCHES 

/INITB/ IBOX(4,6) , IRAY(4,3) ,NRAY,IPT1( 12) , IPT2(12) 
/INPUTA/ EDGE (6, 12) ,POINT(6,8) , SIDE (6, 6) 

/IWRITE/ IWRTC , IWRTI , IWRTN 
/MAPING/ MAPS IDE ( 6 ) ,MAPSEG( 10,12) 

/MARCHS/ MARCH, INDEX( 3) 


/MAXIMUM/ ETAMAX( 10, 12) , SEGMAXt 6 , 10 , 12 ) 

/SPACING/ ISTRTCH( 3) , STRETCH ( 3 ), ETAS ( 3 , 200 ) , ETA( 3 ) , DETA( 3 ) 
/SPLIT/ NSPLIT, ISURF( 10) ,IBC( 10) , 

ISTARTK 10) , INC1 ( 10 ) , ISTART2 ( 10 ) ,INC2(10) 

/UNITS/ NU5,NU6,NU19,NU20,NU21 
/ZONING/ ISECT, NSECT, IZONE,NMBRNDS( 3) 


DIMENSION IFLAGP( 8 ) , IFLAGE( 12) , IFLAGS( 6) 
DIMENSION LSTRTCH ( 3 ) , LBCINPT( 6 ) , LMAPS( 6 ) 
DIMENSION MAPEDGE( 12) ,STRTCH(3) , LNMBRND( 3 ) 


MAPEDGE(I) INDICATES TYPE OF GEOMETRY FOR EDGE I 
= 1 LINEAR 
= 2 CIRCULAR ARC 
= 3 EDGE OF REVOLUTION 


MAPSIDE ( I ) TYPE OF GEOMETRY FOR SURFACE I 

= 1 FLAT SURFACE 
= 4 EDGE OF REVOLUTION 
= 5,6,7 SPECIAL SURFACE SYBROUTINE 

MARCH ETA DIRECTION IN WHICH COMPUTATION IS TO ADVANCE 


NMBRNDS ( I ) NUMBER OF NODES IN EACH ETA(I) DIRECTION 


ISTRTCH (I ) 


0 NO STRETCHING IN ETA(I) DIRECTION 

1 INPUT N VALUES OF ETA(I) 

2 TIGHTEN GRID IN ETA(I) DIRECTION 

3 LOOSEN GRID IN ETA(I) DIRECTION 


WRITE(NU6, 1100) ITITLE, ISECT, NSECT, IZONE 
IF( ISECT. EQ.l) GO TO 100 

SAVE INPUT PARAMETERS IN COMMON WITH PREVIOUS SECTION 


DO 10 1=1,3 

LSTRTCH ( I ) = ISTRTCH ( I ) 
10 LNMBRND(I) = NMBRNDS { I ) 
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C 

DO 20 1=1,6 

20 LMAPS(I) = MAPSIDE(I) 


READ INPUT PARAMETERS 


100 CONTINUE 

READ( NU5 , 10 30 ) MAPEDGE 

READ(NU5, 1000) MAPSIDE 
READ(NU5, 1000) MARCH, NSPLIT, IHGM 

IF( IHGM.EQ.l) THEN 

READ(NU5, 1020) ETA1 1 , ETA12 , ETA13 
READ( NU5, 1020) E : A3 1 , ETA32 , ETA33 
READ( NU5 , 1020 ) STR1 1 , STR1 2 , STR13 
READ(NU5, 1020) STR3 1 , STR32 ,STR3 3 ,STR34 

END IF 

— READ INPUT FOR GENERALIZED BOUNDARY CONDITIONS 
READ( NU5 , 1000 ) NMBRNDS , ISTRTCH 

— PRINT OUT INPUT VARIABLES 

WRITE (NU6,3000) MAPEDGE 
WRITE (NU6, 3010) MAPSIDE 
WRITE (NU6, 3020) MARCH 
WRITE ( NU6 , 3040 ) NMBRNDS , ISTRTCH 

— INITIALIZE FLAGS TO CHECK IF POINT, EDGE, OR SURFACE HAS BEEN INPUT 

DO 130 1=1,8 
130 IFLAGP(I) = 0 

DO 140 1=1,12 
140 IFLAGE(I) = 0 

DO 150 1=1,6 
150 IFLAGS(I) = 0 

— SET ORDER OF EXECUTION 

INDEX(l) = MARCH 
INDEX( 3 ) = 3 

DO 160 1=2 , IDIM 

INDEX(I) = INDEX(I-l) + 1 
160 IF ( INDEX( I ) . GT. IDIM) INDEX(I) = 1 


— INITIALIZE ETA 

ETA(l) =0.0 
ETA( 2 ) =0.0 
ETA( 3 ) = 0.0 
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IF(ISECT.EQ.l ) GO TO 300 

TRANSFER INPUT PARAMETERS IN COMMON WITH PREVIOUS SECTION 

TRANSFER POINTS 

DO 200 1=1, NRAT 

L = IRAT( I, MARCH) 

NEWPT = IPTl(L) 

LASTPT = IPT2( L) 

IFLAGP( NEWPT) = 1 

DO 200 J=1 , 6 

200 POI NT (J, NEWPT) = POINT? J, LASTPT) 

DETERMINE WHICH TWO SURFACES ARE IN COMMON 

NEWSRF = 10 + (MARCH - 6)*MARCH 
LASTSRF = 9 + (MARCH - 7)*MARCH/2 

NUMBER OF EDGES IN COMMON WITH PREVIOUS SECTION 

MEDGES = 3*IDIM - 5 
DO 230 1=1, MEDGES 

DETERMINE WHICH EDGES ARE IN COMMON WITH PREVIOUS SECTION 

NEWEDG = IBOX(I, NEWSRF) 

LASTEDG = IBOXd, LASTSRF) 

I FLAGE( NEWEDG) = 1 

TRANSFER NUMBER OF SEGMENTS 

NMBRSEG( NEWEDG) = NMBRSEG( LASTEDG ) 

TRANSFER SEGMENT MAPPING 

DO 210 J=1,NMBRSEG( LASTEDG) 

MAPSEG(J, NEWEDG) = MAPSEG( J , LASTEDG ) 

TRANSFER EDGE COEFFICIENTS 

DO 210 K= 1 , 8 

COEFE ( K , J , NEWEDG) = COEFE ( K ,J , LASTEDG) 

210 COEFE(K,J, LASTEDG) = 0.0 

I F ( NMBRSEG( LASTEDG) . EQ. 1 ) GO TO 230 

—TRANSFER ETA MAXIMUMS 

DO 220 J=1,NMBRSEG( LASTEDG) - 1 
C 

ETAMAX(J, NEWEDG) = ETAMAX( J , LASTEDG ) 
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ETAMAX( J , LASTEDG) = 0.0 
— TRANSFER SEGMENT MAXIMUMS 
DO 220 K=l,6 

SEGMAX(K,J, NEWEDG) = SEGMAXt K ,J , LASTEDG ) 

220 SEGMAX(K,J, LASTEDG) = 0.0 

—TRANSFER SURFACE MAPPING 

230 MAPSIDE(NEWSRF) = LMAPS ( LASTSRF ) 

IFLAGS(NEWSRF) = 1 

TRANSFER SURFACE COEFFICIENTS 

DO 240 1=1,8 

COEFS ( I , NEWSRF ) = COEFS ( I , LASTSRF ) 

240 COEFS ( I, LASTSRF) = 0.0 

TRANSFER NUMBER OF NODES AND STRECHING DATA FOR ETA_2 AND ETA_3 

DO 250 1=2 , IDIM 

J = INDEX(I) 

NMBRNDS ( J ) = LNMBRND(J) 

STRTCH(J) = STRETCH(J) 

250 ISTRTCH ( J ) = LSTRTCH(J) 


READ INPUT PARAMETERS FOR STRETCHING FUNCTIONS 


300 DO 310 1=1, IDIM 
STRETCH ( I ) = 0.0 
IF( ISTRTCH ( I ) . LE . 1 ) GO TO 310 
READ( NU5 , 1020 ) STRETCH 
GO TO 320 
310 CONTINUE 


COMPUTE STRETCHING FUNCTION PARAMETER B USING NEWTON-RAPHSON 


320 DO 370 1=1, IDIM 

B = STRETCH ( I ) 

DS = B 

TNODE = REAL ( NMBRNDS ( I ) ) 

— DECREASING OR INCREASING SPACING (INPUT MINIMUM SPACING) 

IF( ISTRTCH ( I ) . EQ. 5 .OR. ISTRTCH ( I ) . EQ . 6 ) THEN 
C 

B = SQRT(1.0 - (TNODE - 1 . 0 ) *DS)/( TNODE - 1.0) 
C 

DO 330 ITER=1 , 10 

C - 

ARG1 = B* (TNODE - 1.) 

ARG2 = B* (TNODE - 2.) 
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C 


C 


c 

c 


c 


c 


EXPAl 

= 

EXP(ARGl) 


EXPA2 

= 

EXP{ ARG2 ) 


TANH1 

= 

(EXPAl - 1 ./EXPA1)/(EXPA1 + 

1. /EXPAl) 

TANH 2 

= 

( EXPA2 - 1 ./EXPA2 )/ ( EXPA2 + 

1./EXPA2) 

PSI 

a 

DS - (1.0 - TANH2/TANH1 ) 


TANHP1 

3 

1.0 - TANH1**2 


TANHP2 

=s 

1.0 - TANH2**2 


PSIP 


( 1 . 0/TANH1) * ( TANHP2* { TNODE 

- 2.) 


- ( TANH2/TANH1 ) * TANHP1 *( TNODE - 1.)) 


C 


IF( PSIP . EQ. 0 . ) THEN 


C 

C 

C 


WRITE( NU6 , 5100 ) 

STOP 

END IF 


BO = B 

B = BO - PSI/PSIP 
DBF = (B - B0)/B0 


330 IF(ABS(DBF) . LE . 0.001) GOTO 340 

340 END IF 

DOUBLE STRETCHING (INPUT MINIMUM SPACING) 


C 


IF ( ISTRTCH( I) .EQ.7) THEN 


C 


B = SQRTd.O - ( TNODE - 1 . 0 ) *DS )/ ( TNODE - 1.0) 


DO 350 
C 

ARG1 

ARG3 

C 

EXPA1 

EXPA3 

C 

TANH1 
TANK 3 
C 

PS I 
C 

TANHP1 

TANHP3 

C 

PSIP 

& 

C 


ITER=1 , 10 

= B*( TNODE - 1.) 

= B*( TNODE - 3.) 

= EXP(ARGl) 

= EXP( ARG3 ) 

= ( EXPAl - 1 ./EXPA1 ) / ( EXPAl + 1./EXPA1) 
= ( EXPA3 - 1 ./EXPA3)/(EXPA3 + 1./EXPA3) 

= DS - 0.5* (1.0 - TANH3/TANH1 ) 

= 1.0 - TANH1**2 
= 1.0 - TANH 3**2 

= ( 0 . 5/TANH1 ) * ( TANHP3* ( TNODE - 3.0) 

- (TANH3/TANH1 ) * TANHP1* ( TNODE - 1.0)) 


C 


IF(PSIP.EQ.O.) THEN 

WRITE ( NU6 » 5000 ) 
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C 

C 

C 


STOP 
END IF 

BO = B 

B = BO - PSI/PSIP 
DBF = (B - BO )/B0 


350 IF(ABS(DBF) . LE . 0.001) GO TO 360 

360 END IF 

370 STRETCH ( I ) = B 

• — TRANSFER STRETCHING IN COMMON WITH PREVIOUS SECTION 
IF( ISECT.GT. 1) THEN 

DO 380 1=2 i IDIM 

380 STRETCH ( INDEX( I) ) = STRTCH( INDEX( I) ) 

END IF 

READ INPUT PARAMETERS FOR ARBITRARY GRID SPACING (ETAS) 
400 DO 410 1=1, IDIM 

IF( ISTRTCH ( I) . NE . 1 ) GO TO 410 

IF( ISECT.GT. 1 .AND. I.NE. MARCH) GO TO 410 

READ( NU5 , 1020 ) (ETAS(I,J) , J= 1 , NMBRNDS ( I ) ) 

410 CONTINUE 

READ INPUT PARAMETERS FOR EDGE COEFFICIENTS 

WRITE (NU6, 1100) ITITLE, ISECT,NSECT, IZONE 
WRITE ( NU6 , 2040 ) 

LINE = 6 
II = 0 

■ — TOTAL NUMBER OF EDGES 

NEDGES = 8*IDIM - 12 

DO 540 1=1, NEDGES 

IF( IFLAGE ( I ) . EQ • 1 ) GO TO 540 

ITOTAL = 1 

MAP = MAPEDGE(I) 

— DETERMINE THE NUMBER OF SEGMENTS ON AN EDGE 
DO 500 J=l,10 
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NMBRSEG(I) = J 

MAP = MAP/ 10 
C 

IF(MAP.EQ.O) GO TO 510 
C 

500 ITOTAL = I TOTAL* 10 

C 

510 MAP = MAPEDGE(I) 

C 

DO 530 J=1 , NMBRSEG ( I ) 

C 

DO 520 K=l,8 
520 COEFE ( K , J , I ) = 0.0 

DETERMINE THE MAPPING FOR EACH SEGMENT 

MAPSEG ( J , I ) « MAP/ ITOTAL 

MAP = MAP - MAPSEG( J , I ) * ITOTAL 

— READ IN THE EDGE COEFFICIENTS FOR EACH SEGMENT 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

530 


IF( MAPS EG ( J,I).LE.l) GO TO 530 

READ(NU5, 1020) ( COEFE ( K , J , I ) ,K=1,8) 

IF(I.NE.II) THEN 

LINE = LINE + 2 
I F( LINE .GE . 60 ) THEN 

WRITE(NU6, 1100) ITITLE , I SECT , NSECT , I ZONE 
WRITE ( NU6 , 2040 ) 

LINE = 8 

END IF 

WRITE(NU6, 2050) I , J , ( COEFE ( K , J , I ) ,K=1,8) 

ELSE 

LINE = LINE + 2 
IF( LINE. GE. 60) THEN 

WRITE(NU6, 1100) ITITLE , ISECT , NSECT, IZONE 
WRITE ( NU6 , 2040 ) 

LINE = 8 

END IF 

WRITE(NU6, 2060) J , ( COEFE (K,J,I) ,K=1»8) 

END IF 

II = I 

ITOTAL = I TOTAL/10 
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C 

540 CONTINUE 


READ INPUT PARAMETERS FOR SURFACE COEFFFICIENTS 


LINE = LINE + 5 
IF ( LINE. GE. 60) THEN 

WRITE (NU6, 1100) ITITLE » I SECT , NSECT , I ZONE 
LINE = 6 
END IF 

WRITE ( NU6 , 3060 ) 

DO 610 1=1,6 


C 

c 

c 

c 

c 

c 


IF( I FLAGS ( I) . EQ. 1) GO TO 610 

DO 600 J=1 , 8 
600 COEFS ( J , I ) = 0.0 

I F ( MAPS I DE ( I ) . LE . 2 ) GO TO 610 

READ( NU5 , 1020 ) ( COEFS ( J , I ) , J=1 , 8 ) 

LINE = LINE + 2 

IF( LINE. GE. 60) THEN 

WRITE (NU6, 1100) ITITLE, ISECT, NSECT, IZONE 
WRITE ( NU6 , 3060) 

LINE = 8 

END IF 

WRITE ( NU6, 3065) I, (COEFS( J, I) ,J=1,8) 

610 CONTINUE 

READ INPUT DATA FOR CORNER POINTS AND SEGMENT END POINTS 
LINE = LINE + 5 
IF( LINE.GE .60) THEN 

WRITE (NU6, 1100) ITITLE, ISECT, NSECT, IZONE 
LINE = 6 
END IF 

WRITE ( NU6 , 3070 ) 
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C 


DO 750 1=1, NEDGES 

IF( I .GT. 8 .OR. IFLAGP( I) .EQ. 1) GO TO 700 

READ IN CORNER POINTS 

READ(NU5, 1020) (P0INT{J,I) ,J=1,3) 

LINE = LINE + 2 
IF( LINE .GE . 60 ) THEN 

WRITE (NU6, 1100) I TITLE, ISECT,NSECT, I ZONE 
WRITE ( NU6 , 3070 ) 

LINE = 8 

END IF 

WRITE ( NU6 , 3080 ) I , ( POINT ( J , I ) ,J=1,3) 

700 IF( IFLAGE( I) .EQ. 1 .OR. NMBRSEG ( I ) . EQ . 1 ) GO TO 750 
— READ IN SEGMENT MAXIMUMS 

DO 740 J= 1 , NMBRSEG ( I ) - 1 

READ( NU5 » 1020 ) ( SEGMAX( K , J , I ) , K*1 , 3 ) , ETAMX 
LINE = LINE + 2 
IF( LINE . GE . 60 ) THEN 

WRITE (NU6, 1100) ITITLE, ISECT , NSECT , IZONE 
WRITE ( NU6 , 3070 ) 

LINE = 8 

END IF 


WRITE (NU6, 3090) I , J , ( SEGMAX( K, J, I ) ,K=1,3) , ETAMX 

DO 720 N= 1 , 3 
DO 720 L= 1 , 4 

720 IF{IRAT(L,N) .EQ.I) GO TO 730 

730 CONTINUE 


CONVERT NODE NUMBER TO ETA VALUE 

ETAMAX ( J , I ) = (ETAMX - 1 . 0 ) / ( NMBRNDS ( N ) - 1.0) 


740 CONTINUE 
750 CONTINUE 


RETURN 
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1000 FORMAT ( 1615) 

1010 FORMAT ( 6110) 

1020 FORMAT ( 8G 10 .0 ) 

1030 FORMAT( (4110) ) 

C 

1100 FORMAT ( 1H1,10X,20A4,13X,8H SECTION , 12 , 3H OF,I3,9H FOR ZONE, 13) 
1110 FORMAT (// 11H FIXED ETA_,I1,8H VALUES: // ( 10 ( 3X, F10 . 7 ) ) ) 

2040 FORMAT ( // 25H EDGE SHAPE COEFFICIENTS: 

1 ,// 4 2H. EDGE SEGMENT COEFF_l COEFF_2, 

2 42H COEFF_3 COEFF_4 COEFF_5, 

3 42H COEFF_6 C0EFF_7 COEFF_8 ) 

2050 FORMAT (/ 2X, 12 , 6X, II , 3X, 8 ( IX, FI 3 . 7) ) 

2060 FORMAT ( / 10X, 1 1 , 3X, 8 ( IX, FI 3 . 7 ) ) 

2080 FORMAT ( / 7H POINT , 1 2 , 1H : , 1 2X, 2 ( 7X, F13 . 7 ) , 9X, F7 . 2 ) 

2090 FORMAT ( / 7H EDGE , 12 , 1H : , 7X, 13 , 2X, 2 ( 7X, F13 . 7 ) , 9X, F7 . 2 , 9X, F9 . 2) 


3000 


3010 


EDGE_4, 
EDGE_8 , 
EDGE 12 


FORMAT (// 32H EDGE SHAPE FUNCTION INDICATORS: 

// 40H EDGE_1 EDGE_2 EDGE_3 

40H EDGE_5 £DGE_6 EDGE_7 

40H EDGE_9 EDGE_10 EDGE_1 1 

FORMAT!// 35H SURFACE SHAPE FUNCTION INDICATORS: 

// 45H SURFACE_1 SURFACE_2 SURFACE 

45H SURFACE_4 SURFACE_5 SURFACE" 

/ 6( 10X, 15) ) 

FORMAT!/// 26H MARCHING DIRECTION = ETA_,I1) 

FORMAT!// 17H NUMBER OF NODES: 

30H ETA_1 ETA_2 ETA_3 / 3110 

29H NODAL STRETCHING INDICATORS: 

30H ETA_1 ETA_2 ETA_3 / 3110 ) 

32H STRETCHING FUNCTION PARAMETERS: 

45H ETA_1 ETA_2 

3 ( 2X f F13 . 7) ) 

28H SURFACE SHAPE COEFFICIENTS: 

42H SURFACE COEFF_l COEFF 

42H COEFF_3 COEFF_4 COEFF* 

4 2H COEFF_6 COEFF_7 COEFF*_ 

FORMAT! / 4X,I1,9X,8( 1X.F13.7) ) 

FORMAT (//' COORDINATES 

//10X,' SEGMENT X T' , 

' Z ETAMAX ' ) 

FORMAT!/ 7H POINT , 12, 1H: , 8X, 3 ( 6X, F13 . 7) ) 

FORMAT!/ 7H EDGE , 12 , 1H: , 3X, 13 , 2X, 3 ( 6X, F13 .7 ) 


/ 12110 ) 


3, 

'6 


3020 
3040 

1 // 
2 // 
3 // 

3050 FORMAT!// 
1 // 
2 / 
3060 FORMAT!// 
1 // 
2 
3 

3065 
3070 

1 
2 

3080 
3090 

8 

5000 
5100 


ETA 3 


2 , 
"5, 
*8 ) 


FORMAT! 1H0, 
FORMAT! 1 HO, 1 


DEAD 

DEAD 


IN 

IN 


INPUT* 

INPUT' 


) 


,7X,F9.2) 


END 




SUBROUTINE GRID 


GRID CONTROLS THE GENERATION OF THE SPATIAL COORDINATES OF EACH NODE 
IN THE GRID. 


C 


DOUBLE PRECISION NODE 
REAL NODE 
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C 


C 


C 


COMMON /COUNTER/ 
COMMON /HEADER/ 
COMMON /INITA/ 
COMMON /INITB/ 
COMMON /INPUTA/ 
COMMON /INPUTBC/ 
COMMON /IWRITE/ 
COMMON /MAPING/ 
COMMON /MARCHS/ 
COMMON /NODNMBR/ 
COMMON /OUT/ 
COMMON /PARTIAL/ 
COMMON /SPACING/ 
COMMON /UNITS/ 
COMMON /ZONING/ 


NODESAV , NODETOT , NBNODES , NPLANE 
ITITLE ( 20 ) /LINE 
IDIM,MAPTEN, INCHES 

IBOX(4,6) ,IRAY(4,3) ,NRAY, IPT1 ( 12) , IPT2( 12) 

EDGE( 6/12) /POINT (6»8) /SIDE (6,6) 

ISIDE( 3) 

IWRTC r IWRTI , IWRTN 
MAPS IDE ( 6 ) ,MAPSEG{ 10,12) 

MARCH, INDEX( 3) 

NODENUM( 10000) ,TOL( 3) 

NODE( 5, 10000) 

DEDN( 3,12) , DSDN( 3,2) , SNORMAL( 3,6) 

ISTRTCH( 3 ) , STRETCH ( 3 ) , ETAS (3,200) ,ETA( 3) ,DETA( 3) 
NU5 , NU6 ,NU19 , NU20 , NU21 
I SECT , NSECT , I ZONE , NMBRNDS ( 3 ) 


DIMENSION E ( 26 ) 

DIMENSION I FACE 1( 3,3) , IFACE2(3,3) ,ISIDES(6) 
DIMENSION TANGENT ( 3) ,VECTOR( 3) ,DIRECT(3) 


DATA IFACE1 /0 , 1 , 2 , 1 , 0 , 5 , 2 , 5 , 0/ 
DATA IFACE2 /0, 3, 4, 3, 0,6, 4, 6,0/ 
DATA ISIDES /3 , 4 , 1 , 2, 6 , 5/ 


INIT = 1 


DO 30 1=1, IDIM 

IF(STRETCH(I) .LT.1.0) STRETCH ( I ) = 1.0 

DETA(I) = 1 . 0/ ( NMBRNDS ( I ) - 1.0) 

FACTOR = NMBRNDS ( I) *10.0 
IF ( ISTRTCH ( I ) . NE . 0 ) FACTOR = NMBRNDS ( I )* 20 . 0 

DO 20 J=1,NRAY 

L = I RAT ( J , I ) 

• — COMPARE X, Y , AND Z BETWEEN END POINTS OF AN EDGE 
DO 20 K= 1 , IDIM 

DELTA = ABS ( POINT ( K , IPT1 ( L ) ) - POINT( K , IPT2 ( L) ) )/FACTOR 

IF( DELTA. LE . 0 . 0 ) GO TO 20 

IF(DELTA.LT.TOL(K) ) TOL(K) = DELTA 

20 CONTINUE 
30 CONTINUE 

TOL(l) = .005 
TOL( 2 ) = .005 
TOL( 3 ) = .005 

— DETERMINE SURFACES 

ISIDE1 = IFACE 1 ( MARCH, INDEX( 2 ) ) 

ISIDE2 = I FACE 2 ( MARCH , INDEX( 2 ) ) 

ISIDE3 = IFACE 1 ( MARCH , INDEX( 3 ) ) 
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) 

ISIDE4 = I FACE 2 ( MARCH , INDEX( 3 ) ) 

ISIDE5 = IFACE1 ( INDEX( 2 ) , INDEX ( 3 ) ) 

ISIDE6 = I FACE 2 ( INDEX ( 2 ) , INDEX ( 3 ) ) 

MAXIMUM NUMBER OF NODES IN A PLANE 

MAXPL = NMBRNDS( INDEX( 2) ) *NMBRNDS( INDEX( 3) ) 

INDEX(l) = MARCH 

AXIS 

DO 700 I AXIS- 1 , NMBRNDS ( INDEX( 1) ) 

I F ( I AXIS . EQ . 1 .AND. ISECT.GT.l) GO TO 700 

SEPERATE BOUNDARY CONDITIONS AND DETERMINE ETA 

CALL ETABC ( MARCH , INDEX ( 1 ) , I AXIS) 

CALCULATE COORDINATES AND DERIVATIVES FOR POINTS ON EDGES 

CALL EDGES( INIT, INDEX( 1 ) , ETA( INDEX( 1 ) ) ) 

NUMBER OF NODES TO BE STORED 

NODSTOR = NODESAV 

ROW 

DO 500 J AXIS 3 1 , NMBRNDS (INDEX(2)) 

SEPERATE BOUNDARY CONDITIONS AND DETERMINE ETA 

CALL ETABC(MARCH, INDEX( 2) ,JAXIS) 

CALCULATE COORDINATES AND DERIVATIVES FOR POINTS ON EDGES 

CALL EDGES ( INIT, INDBX( 2 ) , ETA( I NDEX( 2 ) ) ) 

— CALCULATE COORDINATES AND SURFACE NORMAL FOR POINTS ON SURFACES 
CALL SURFACE (INIT,ISIDE1) 

CALL SURFACE( INIT, ISIDE2) 

COLUMN 

DO 400 KAXIS=1 , NMBRNDS ( INDEX( 3 ) } 

— SEPERATE BOUNDARY CONDITIONS AND DETERMINE ETA 
CALL ETABC ( MARCH , INDEX ( 3 ) , KAXIS ) 

— CALCULATE COORDINATES AND DERIVATIVES FOR POINTS ON EDGES 
CALL EDGES( INIT, INDEX( 3) ,ETA( INDEX( 3) ) ) 
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C CALCULATE COORDINATES AND SURFACE NORMAL FOR POINTS ON SURFACES 

C 

CALL SURFACE ( INIT,ISIDE3) 

C 

CALL SURFACE( INIT, IS IDE 4) 

C 

CALL SURFACE( INIT, ISIDE5) 

C 

CALL SURFACE( INIT, ISIDE6) 

C 

C ETA COEFFICIENTS FOR TRI-LINEAR INTERPOLATION 


Ed) 


1.0 - ETA( 3 ) 

E(2) 

= 

ETA( 3) 

E(3) 

= 

1.0 - ETA( 2 > 

E ( 4 ) 

= 

ETA( 2) 

E( 5) 

S 

1.0- ETA{1) 

E ( 6 ) 

= 

ETA( 1) 

E(7) 

s 

E( 5)*E( 3) 

E ( 8 ) 

= 

E ( 5 ) *ETA( 2) 

E ( 9 ) 

a 

ETA( 1 ) *E( 3 ) 

E ( 10) 

= 

ETA( 1 ) * ETA( 2 ) 

E ( 11) 

= 

E ( 5 ) *E( 1 ) 

E( 12) 

= 

E ( 5 ) *ETA( 3 ) 

E( 13 ) 

a 

ETA( 1 ) * E( 1 ) 

E( 14) 

3 

ETA( 1 ) * ETA( 3 ) 

E ( 15) 

= 

E ( 3 ) * E ( 1 ) 

E ( 16 ) 

= 

E ( 3 ) *ETA( 3 ) 

E ( 1 7 ) 


ETA( 2 ) * E ( 1 ) 

E( 18 ) 


ETA( 2 ) *ETA( 3) 

E( 19) 

= 

E( 5 ) * E( 3 ) *E( 1 ) 

E ( 20 ) 

S 

E( 5 ) *E( 3 ) *ETA( 3 ) 

E( 21 ) 

= 

E ( 5 ) * ETA ( 2 ) * E ( 1 ) 

E ( 22 ) 

= 

E(5)*ETA(2)*ETA(3) 

E( 23 ) 

= 

ETA( 1 ) * E ( 3 ) *E ( 1 ) 

E( 24 ) 

= 

ETA( 1 ) ‘ E ( 3 ) *ETA( 3 ) 

E( 25 ) 

= 

ETA ( 1)*ETA( 2)*E( 1) 

E ( 26 ) 

= 

ETA( 1 )*ETA( 2) *ETA( 3 ) 


C INCREMENT NODE COUNTERS 

C 

INIT = 0 


C 


NODESAV = NODESAV + 


1 


» 


C 


NODNUM = NODESAV + NODETOT 


C CALCULATE THE COORDINATES 

C 

DO 340 L=1 , 3 

340 NODE(L, NODESAV) = E( 1 ) *SIDE{ L, 1 ) + E( 2)*SIDE( L, 3) 

1 + E ( 3 ) ‘SIDE ( L, 2 ) + E(4)*SIDE(L,4) 

2 + E( 5) *SIDE( L, 5 ) + E( 6 ) ‘SIDE ( L , 6 ) 

3- E( 7 ) *EDGE( L, 5 ) -E ( 8 ) *EDGE(L, 8)-E{ 9 ) *EDGE ( L, 6 ) -E ( 10 ) ‘EDGE ( L, 7) 

4- E( 11)*EDGE(L,4)-EU2)*EDGE(L,12)-E(13)*EDGE(L,2)-E( 14 ) *EDGE( L, 10 ) 

5- E( 15 ) * EDGE ( L»l)-E(16) ‘EDGE ( L, 9)-E( 17 ) ‘EDGE ( L, 3 ) -E( 18)*EDGE(L, 11) 

6 + E( 19 ) ‘POINT! L , 1 ) + E ( 20 ) * POI NT ( L , 5 ) 

7 + E( 21)* POI NT (L, 4) + E( 22)* POI NT (L, 8) 

8 + E(23)*POINT(L,2) + E( 24 ) * POINT( L, 6 ) 

9 -t- E ( 25 ) * POINT ( L, 3 ) + E ( 26 ) *POINT{ L, 7 ) 

C 
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C- 

C 

C 

C 

C 

C 


— INTERMEDIATE ERROR PRINT 

IF( IWRTI.EQ.O .OR. MOD( NODNUM, IWRTI ) . NE.O) GO TO 350 
LINE = LINE + 27 
IF( LINE. GE. 60) THEN 

WRITE (NU6, 1100) ITITLE, ISECT, NSECT, I ZONE 
LINE = 28 


C 

C 


C 


C 


C 


C 

C 

C 


END IF 


WRITE(NU6,1110) NODNUM 

1 ,E( 1) , ( SIDE( L, 1 ) , L=1 , 3) 

2 ,E(3) ,(SIDE(L,2) ,L=1,3) 

3 ,E(2),(SIDE(L,3),L=1,3) 

4 ,E(4) , ( SIDE ( L , 4 ) , L=1 , 3 ) 

5 ,E(5) , ( S I DE ( L , 5 ) # L= 1 # 3 ) 

6 ,E(6) , (SIDE(L,6) ,L=1 ,3) 


WRITE (NU6, 1120) E(15) .(EDGE (L,l) , L= 1,3) 

1 / E( 13 ) , ( EDGE ( L, 2 ) < L= 1 , 3 ) 

2 , E( 17 ) , ( EDGE ( L, 3 ) » L*1 * 3 ) 

3 , E( 1 1 ) , ( EDGE ( L* 4 ) , L=1 , 3) 

4 / E ( 7) r (EDGE(L,5) ,L=1,3) 

5 ,E( 9 ) , ( EDGE ( L , 6 ) f L a 1 , 3 ) 


WRITE (NU6, 1130) E( 10 ) , ( EDGE ( L , 7 ) , L= 1 , 3 ) 

1 ,E( 8 ) , ( EDGE ( L , 8 ) , L= 1 , 3 ) 

2 , E( 16 ) , ( EDGE ( L» 9 ) , L= 1 , 3 ) 

3 , E ( 14 ) , ( EDGE (L r 10) f L=l,3) 

4 , E( 18 ) , ( EDGE (L,ll) ,L=1,3) 

5 ,E{ 12) , (EDGE(L,12) ,L=1,3) 


1 

2 

3 

4 

5 

6 
7 


WRITE (NU6, 1140) E ( 19 ) , ( POINT( L , 1.) ,L=1,3) 
, E ( 23 ) ,( POINT (L,2) ,L=1,3) 
f E( 25 ) i ( POINT ( L» 3 ) , L»1 » 3) 
, E( 21 ) , ( POINT( L, 4 ) ,L=1,3) 
, E( 20 ) , { POINT( L , 5) , L=1 * 3 ) 
,E( 24 ) , ( POINT ( L, 6) , L= 1 » 3 ) 
, E( 26 ) , ( POINT{ L , 7) , L=1 » 3 ) 
, E( 22) , ( POINT ( L, 8 ) , L=1 , 3 ) 


350 CONTINUE 


NODE( 4 , NODESAV) = THETA 
NODE ( 5 , NODESAV ) = PHI 
C 

400 CONTINUE 
500 CONTINUE 


C PLANE OUTPUT 


IF(IAXIS.EQ.l) GO TO 700 
I F ( NODSTOR .NE.O) THEN 
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C 

C 

C 

C 

C 

C 

C 

C- 

C 

C- 

C 

c- 

c 

c 

c 

c 


■ — PRINT AND STORE DATA 

CALL OUTPUT ( NU20 , NODSTOR) 

• — TRANSFER DATA SECOND PLANE TO FIRST PLANE 
DO 630 1=1 , NODESAV 
DO 600 J=1 , 5 

600 NODE ( J , I ) = NODE ( J / I + NODSTOR) 

630 CONTINUE 
END IF 

—TRANSFER NODE NUMBERS FROM SECOND PLANE TO FIRST PLANE 
DO 650 L=1,MAXPL 

650 NODENUM(L) = NODENUM( L + MAXPL) 

700 CONTINUE 
SECTION OUTPUT 

IF( ISECT.NE.NSECT) RETURN 
■ — PRINT AND STORE DATA 


NODSTOR = NODESAV 

IF( NODESAV. NE.0) CALL OUTPUT ( NU20 , NODSTOR) 

NODESAV = 0 


RETURN 


C 

C FORMAT STATEMENTS 

C 


1100 FORMAT ( 1H1, 10X,20A4, 13X,8H SECTION , 12 , 3H OF,I3,9H FOR ZONE, 13) 
C 


1110 FORMAT ( 


8H 

NODE 

= ,110 




1 

/ 

41X,2H X 

, 13X, 2H Y,12X,2H 

Z , 20X 

, 1 5H 

1 FLOW DIRECTION 

2 

/ 

8H 

E( 1) 

=,F13.7, 15H 

SIDE(l) 

=,6(1X,F13.7) 

3 

/ 

8H 

E( 3) 

=,F13.7 , 15H 

SIDE ( 2 ) 

= , 6 ( IX, F13 . 7 ) 

4 

/ 

8H 

E( 2) 

* ,F13 . 7 , 15H 

SIDE( 3) 

= , 6 ( IX, F13 . 7 ) 

5 

/ 

8H 

E< 4) 

=,F13.7, 15H 

SIDE( 4 ) 

=,6(1X,F13.7) 

6 

/ 

8H 

E( 5) 

=,F13.7,15H 

SIDE( 5) 

=,6(1X,F13.7) 

7 

/ 

8H 

E( 6) 

=,F13.7,15H 

SIDE ( 6 ) 

=,6(1X,F13.7) ) 

1120 FORMAT ( 


8H 

E ( 1 5 ) 

=,F13.7,15H 

EDGE ( 

1) 

=,6(1X,F13.7) 

1 

/ 

8H 

E( 13) 

=,F13.7, 15H 

EDGE ( 

2) 

= ,6( IX, FI 3 . 7 ) 

2 

/ 

8H 

E( 17) 

=,F13 . 7, 15H 

EDGE ( 

3) 

=,6( 1X,F13.7) 

3 

/ 

8H 

E( 11) 

= , FI 3 . 7 , 15H 

EDGE ( 

4) 

= , 6 ( IX, FI 3 . 7 ) 

4 

/ 

8H 

E { 7) 

=,F13.7, 15H 

EDGE ( 

5) 

=,6(1X,F13.7) 

5 

/ 

8H 

E( 9) 

= , F13 . 7 , 15H 

EDGE( 

6) 

= , 6 ( IX, F13 . 7 ) ) 

1130 FORMAT ( 


8H 

E ( 10 ) 

= ,F13 .7 , 15H 

EDGE( 

7) 

= , 6 { IX, F 13 . 7 ) 

1 

/ 

8H 

E( 8) 

=,F13.7, 15H 

EDGE ( 

8) 

=,6( 1X,F13.7) 

2 

/ 

8H 

E( 16) 

=,F13.7,15H 

EDGE( 

9) 

=,6( 1X,F13.7) 
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3 

/ 

8H 

E( 14 ) 

= , FI 3 . 7 , 15H 

EDGE( 10) 

=,6( 1X,F13.7) 

4 

/ 

8H 

E( 18 ) 

= ,F13 .7, 15H 

EDGE( 11) 

=, 6( IX, F13 . 7 ) 

5 

c 

/ 

8H 

E ( 12 ) 

= , F13 . 7 , 15H 

EDGE( 12) 

= , 6{ IX, FI 3 . 7 ) ) 

1140 FORMAT ( 


8H 

E ( 19 ) 

= , FI 3 . 7 , 15H 

POINT( 1) 

= , 6 ( 1X,F13 . 7 ) 

1 

/ 

8H 

E ( 23 ) 

= , FI 3 . 7 , 15H 

POINT ( 2 ) 

= , 6{ IX, F13 .7) 

2 

/ 

8H 

E ( 25 ) 

= , F13 . 7 , 15H 

POINT ( 3 ) 

=,6( 1X,F13. 7) 

3 

/ 

8H 

E ( 2 1 ) 

= , FI 3 . 7 , 15H 

POINT( 4 ) 

= , 6 ( IX, F13 . 7 ) 

4 

/ 

8H 

E( 20 ) 

= , FI 3 . 7 , 15H 

POINT ( 5 ) 

= , 6 ( 1X,F13 . 7 ) 

5 

/ 

8H 

E ( 24 ) 

= , FI 3 . 7 , 15H 

POINT( 6 ) 

= , 6 ( IX, F13 . 7 ) 

6 

/ 

8H 

E ( 26 ) 

= , F13 . 7 , 15H 

POINTf 7) 

=,6(1X,F13.7) 

7 

/ 

8H 

E( 22 ) 

= , FI 3 . 7 , 15H 

POINT( 8 ) 

=,6(1X,F13.7) ) 


C 


END 


SUBROUTINE EDGES ( INIT, IDIR, ETA) 


INTERPOLATES ALONG EACH EDGE OF EACH SIDE OF ZONE AS PART OF THE 
BI/TRI- LI NEAR INTERPOLATION SCHEME 


COMMON /COEFF/ COEFE ( 8 , 1 0 , 1 2 ) , COEFS ( 8 , 6 ) , NMBRSEG ( 1 2 ) 

COMMON /COUNTER/ NODESA V , NODETOT , NBNODES , NPLANE 
COMMON /HEADER/ ITITLE ( 20 ) , LINE 
COMMON /INITA/ IDIM, MAPTEN, INCHES 

COMMON /INITB/ IBOX( 4 , 6 ) , IRAY ( 4 , 3 ) , NRAY , IPT1 ( 1 2 ) , IPT2 ( 1 2 ) 

COMMON /INPUTA/ EDGE ( 6 , 1 2 ) , POINT{ 6 , 8 ) , SIDE ( 6 , 6 ) 

COMMON /IWRITE/ IWRTC , IWRTI , IWRTN 

COMMON /MAPED/ KSEG< 1 2 ) , UAX IS ( 3 , 6 ) , IEDGE1 { 6 ) , I EDGE2 ( 6 ) , GAMMA 
COMMON /MAXIMUM/ ETAMAX( 10,12) , SEGMAX( 6,10,12) 

COMMON /PARTIAL/ DEDN ( 3 » 1 2 ) , DSDN ( 3 , 2 ) , SNORMAL (3,6) 

COMMON /SEGMENT/ PT1 ( 6 , 10 , 1 2 ) , PT2 ( 6 , 10 , 1 2 ) , ETAl (10,12), ETA2( 10,12) 
COMMON /UNITS/ NU5 , NU6 , NU19 , NU20 , NU2 1 

COMMON /ZONING/ ISECT , NSECT, IZONE , NMBRNDS ( 3 ) 


INTERMEDIATE PRINT 


NODNUM = NODESAV + NODETOT -t- 1 

I F ( IWRTI. EQ.O .OR. MOD ( NODNUM , IWRTI ). NE . 0 ) GO TO 20 
LINE = LINE + 2 
IF( LINE. GE. 60) THEN 

WRITE { NU6 , 1 1 00 ) ITITLE , ISECT , NSECT , IZONE 
LINE = 3 


END IF 


CALCULATE EDGE COORDINATES AND DERIVATIVE 


20 DO 200 1=1, NRAY 
C 
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— DETERMINE WHICH EDGE 

IEDGE = I RAY ( I , IDIR) 

EDGE INITIALIZATION 

IF( INIT.EQ. 1) CALL EDGMAP1 ( IEDGE , I DI R , I ) 

— DETERMINE WHICH SEGMENT 

DO 100 JSEG=1 ,NMBRSEG( IEDGE) 

ISEG = JSEG 

100 IF(ETA.GE.ETA1( JSEG, IEDGE) .AND. ETA. LE. ETA2 ( JSEG , IEDGE ) ) GOTO 110 

110 KSEG( IEDGE) = ISEG 

— DETERMINE WHERE ALONG THE SEGMENT 

DENOM = ETA2( ISEG, IEDGE) - ETA1 ( ISEG , IEDGE ) 

I P ( DENOM.EQ.0.0) DENOM = 1.0 

RATIO = (ETA - ETA1{ ISEG, IEDGE) )/DENOM 

—CALCULATE THE COORDINATES AND DERIVATIVE 

CALL EDGMAP2( IEDGE , ISEG , RATIO , EDGE ( 1 , IEDGE) , DEDN ( 1 , IEDGE) ) 

INTERMEDIATE PRINT 

IF( IWRTI .EQ.O .OR. MOD( NODNUM , IWRTI ) . NE . 0 ) GO TO 200 
LINE = LINE + 1 
IF( LINE. GE. 60) THEN 

WRITE (NU6, 1100) ITITLE , I SECT, NSECT , I ZONE 
LINE = 2 
END IF 

WRITE (NU6, 1120) NODNUM, IEDGE, 

1 (EDGE(J, IEDGE) ,J=1,3) , THETA, PHI, 

2 ETA, RATIO, ETAMAX( ISEG, IEDGE) 

200 CONTINUE 
RETURN 

— FORMAT STATEMENTS 

1100 FORMAT ( 1H1, 10X,20A4,13X,8H SECTION, 1 2 , 3H OF,I3,9H FOR ZONE, 13) 

1120 FORMAT ( IX, 16 , 4X, 12 , 3X, 3 ( 3X, FI 3 . 7 ) , 1 IX, 2 ( 2X, F7 . 2 ) , 3 ( 2X, F7 . 5 ) ) 

END 

t********************************************************************** 
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q * * * **MAPP ING* ** * ** ** * ************************************************** 
£*********************************************************************** 

C 

SUBROUTINE EDGMAP1 ( IEDGE , IDIR , NMBREDG ) 

EDGE MAPPING INITIALIZATION 


COMMON 

COMMON 

COMMON 

& 

COMMON 

& 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


/COEFF/ COEFE ( 8,10,12), COEFS (8,6), NMBRSEG( 1 2 ) 

/COUNTER/ NODESAV , NODETOT , NBNODES , NPLANE 
/EDGEO/ UI( 3, 10, 12) ,UJ(3,10, 12) , UK (3, 10, 12) , 
R1(3,10,12),R2(3,10,12) , THETA( 10,12) 

/EDGE3/ ARC( 10,12) ,ARC1( 10,12) ,XLENGTH( 10, 12) , 

RA( 10, 12) ,RC( 10, 12) ,RE( 10,12) ,THETA1 ( 10,12) 
/EDGE8/ RM1( 10, 12) ,RM2( 10 , 1 2 ) , RPA1 ( 10 , 1 2 ) ,RPA2(10, 12) 
/HEADER/ ITITLE( 20 ) , LINE 
/INITA/ IDIM,MAPTEN, INCHES 

/INITB/ IBOX( 4,6) , IRAK 4 , 3 ) , NRAT , IPT1 ( 12 ) ,IPT2(12) 

/INITC/ PI , RADDEG 

/INPUTA/ EDGE (6,12) ,POINT<6,8) ,SIDE{6,6) 

/MAPED/ KSEG ( 1 2 ) , UAXIS( 3,6),IEDGE1(6), I EDGE 2 ( 6 ) , GAMMA 
/MAPING/ MAPSIDE ( 6 ) ,MAPSEG( 10, 12 ) 

/MARCHS/ MARCH, INDEX( 3) 

/MAXIMUM/ ETAMAX( 10, 12) , SEGMAX( 6,10,12) 

/PARTIAL/ DEDN ( 3 , 1 2 ) , DSDN ( 3 , 2 ) , SNORMAL (3,6) 

/SEGMENT/ PT 1(6, 10, 12), PT2 (6,10,12) ,BTA1( 10,12) ,ETA2( 10, 12) 


DIMENSION ETAMX( 10,12) ,VECTER( 3) 


INITIALIZE 

KSEG( IEDGE) = 1 
DO 10 N=1,NMBRSEG( IEDGE) 

10 ETAMX(N, IEDGE) = ETAMAX( N , IEDGE ) 


ETAMX(NMBRSEG( IEDGE) , IEDGE) =1.0 

ETA1( 1, IEDGE) = 0.0 
ETAMAX(NMBRSEG( IEDGE) , IEDGE) = 1.0 

INITIALIZE SEGMENTS 


DO 1000 ISEG=1,NMBRSEG( IEDGE) 


— DETERMINE THE COORDINATES AT THE END POINTS 


IF(ISEG.EQ.l) THEN 
C 

DO 50 J= 1 , 6 

PT1(J, 1, IEDGE) 
SEGMAX( J ,NMBRSEG( IEDGE) , IEDGE) 
50 PT2(J,ISEG, IEDGE) 

C 


DO 60 J= 1 , 6 

PT1(J,ISEG, IEDGE) 
60 PT2(J,ISEG, IEDGE) 

C 


POINT( J , IPT1 ( IEDGE) ) 
POINT(J,IPT2( IEDGE)) 
SEGMAX( J , ISEG, IEDGE) 


SEGMAX( J, ( ISEG-1 ) , IEDGE) 
SEGMAX( J, ISEG, IEDGE) 
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END IF 

— ETA VALUE AT THE END OF EACH SEGMENT 

IF(ISEG.GT.l) ETA I { ISEG, I EDGE) = ETA2( ISEG-1 , I EDGE) 

ETA2( ISEG, IEDGE) = ETAMX( ISEG, IEDGE) 

VECTORS FROM ORIGIN TO END POINTS 

DO 80 J= 1,3 

R1 ( J, ISEG, IEDGE) = PT1 (J , ISEG, IEDGE ) - COEFE(J, ISEG, IEDGE) 

80 R2(J, ISEG, IEDGE) = PT2(J, ISEG, IEDGE) - COEFE(J, ISEG, IEDGE) 

CALL VMAG(R1(1, ISEG, IEDGE) ,RMAG1) 

CALL VMAG(R2(1, ISEG, IEDGE) ,RMAG2) 

CHOOSE MAPPING FUNCTION 

MAP = MAPSEGf ISEG, IEDGE) 

GO TO (1000,200,300,400,500) MAP 

CIRCULAR ARC MAP = 2 

SWEEP ANGLE 

200 CALL VDOT(Rl(l, ISEG, IEDGE), R2(l, ISEG, IEDGE) ,RlDOTR2) 

ARG = R 1 DOTR2/ ( RMAG 1 * RM AG 2 ) 

I F ( ABS ( ARG ) . GT .1.0) ARG = ARG/ABS ( ARG) 

THETA( ISEG, IEDGE) = ACOS(ARG) 

GO TO 1000 

EDGE OF REVOLUTION MAP = 4 

300 CONTINUE 

UK: NORMALIZED VECTOR ALONG AXIS OF REVOLUTION 

CALL VMAG ( COEFE ( 4 , ISEG , IEDGE ) , UKM ) 

DO 410 1=1,3 

410 UK(I, ISEG, IEDGE) = COEFE(I + 3 , ISEG , IEDGE )/UKM 

UJ: NORMALIZED VECTOR FROM AXIS TOWARD EDGE 

CALL CROSS ( UK ( 1, ISEG, I EDGE) , R1 ( 1 , ISEG, IEDGE) , UJ ( 1 , ISEG , IEDGE ) ,11) 

UI: NORMALIZED VECTOR FROM AXIS TO FIRST POINT 

CALL CROSS (UJ(1, ISEG, IEDGE) ,UK( 1, ISEG, IEDGE) , UI ( 1 , ISEG , IEDGE ) ,12) 

— Rl: VECTOR FROM AXIS TO FIRST POINT 

CALL VDOT(Rl( 1, ISEG, IEDGE ), UK ( 1, ISEG, IEDGE ) ,RPA1( ISEG, IEDGE) ) 
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CALL VADD( 1 . 0 , R1 ( 1/ IS EG, IEDGE) , -RPA1 ( IS EG/ I EDGE) , UK ( 1 , ISEG , IEDGE) , 
& R1{1, ISEG, IEDGE) / VECTER) 

C 

CALL VMAG(R1(1, ISEG, IEDGE) ,RMl ( ISEG, IEDGE) ) 

R2 : VECTOR FROM AXIS TO SECOND POINT 

CALL VDOT(R2(l, ISEG, IEDGE) ,UK< 1 , ISEG, IEDGE) ,RPA2(ISEG, IEDGE) ) 

CALL VADD(1.0,R2( 1, ISEG, IEDGE ) ,-RPA2( ISEG, IEDGE) , UK ( 1 , ISEG , IEDGE ) , 
& R2(l, ISEG, IEDGE) , VECTER) 

CALL VMAG ( R2 (1, ISEG, IEDGE ) ,RM2( ISEG, I EDGE) ) 

SWEEP ANGLE 

CALL VDOT( Rl ( 1 , ISEG, IEDGE) ,R2( 1 , I SEG , IEDGE ) , R1DOTR2 ) 

ARG = RIDOTR2/ ( RM1 { ISEG , IEDGE )*RM2(IS EG, IEDGE) ) 

IF(ABS(ARG) .GT.l. ) ARG = ARG/ABS( ARG) 

THETA ( ISEG, IEDGE) = ACOS(ARG) 

GO TO 1000 

400 CONTINUE 


USER SUPPLIED SPECIAL EDGE INITIALIZATION MAP = 8 


GO TO 1000 


USER SUPPLIED SPECIAL EDGE INITIALIZATION MAP = 9 


500 CONTINUE 

1000 CONTINUE 

RETURN 

END 

SUBROUTINE EDGMAP2 ( IEDGE , ISEG , RATIO, POINT, TANGENT) 


THIS ROUTINE CALCULATES THE COORDINATES 
OF A POINT ON AN EDGE. 


COMMON /COEFF/ 
COMMON /EDGE0/ 

& 

COMMON /EDGE3/ 

& 


COEFE (8,10,12) , COEFS (8,6) , NMBRSEG( 1 2 ) 

UI(3, 10,12) ,UJ(3,10,12) , UK (3, 10, 12) , 
R1(3,10,12),R2(3,10,12) , THETA { 10,12) 

ARC (10,12) ,ARC1( 10, 12) ,XLENGTH( 10,12) , 
RA(10,12) ,RC( 10, 12) ,RE(10,12) ,THETA1( 10,12) 
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C 


C 

C 


COMMON /EDGE8/ RM1 ( 10 , 12) , RM2( 10 , 12 ) ,RPA1 ( 10 , 12 ) , RPA2( 10 , 12) 
COMMON /MAPING/ MAPSIDE ( 6 ) ,MAPSEG ( 1 0 , 1 2 ) 

COMMON /SEGMENT/ PT1 ( 6 , 10 , 12 ) , PT2 ( 6 , 10 , 12 ) , ETA1 ( 10 , 12 ) , ETA2 ( 10 , 1 2 ) 
COMMON /ZONING/ ISECT , NSECT , IZONE , NMBRNDS ( 3 ) 

DIMENSION POINT( 6 ) , TANGENT( 3 ) 

DIMENSION FLOWI ( 3 ) ,FLOWJ(3) ,FLOWK(3) 

DIMENSION VECTOR* 3) ,C( 3) 

MAP = MAPSEG( ISEG, IEDGE) 

GO TO (100,200,300,400,500) MAP 


LINEAR EDGE 


MAP = 1 


100 DO 110 1=1,3 

POINT(I) = (1. - RATIO)*Rl(I, ISEG, IEDGE) + RATIO* R2 ( I , ISEG , IEDGE ) 
110 VECTOR(I) = R2( I , ISEG, IEDGE) - R1 ( I , ISEG, IEDGE) 

CALL VMAG ( VECTOR , RMAG ) 

IF( RMAG. EQ. 0.0) RMAG = 1.0 
DO 120 1=1,3 

120 TANGENT ( I ) = ( R2 ( I , ISEG , IEDGE ) - R1 ( I , ISEG , IEDGE )) /RMAG 
GO TO 1000 


CIRCULAR ARC MAP = 2 


200 PHI1 = (1.0 - RATIO) *THETA( ISEG, IEDGE) 

PHI 2 = RATIO*THETA( ISEG, IEDGE) 

DO 210 1=1,3 

POINT(I) = COEFE(I, ISEG, IEDGE) + 

1 (SIN(PHI1)*R1(I, ISEG, IEDGE) + SIN ( PHI 2 ) *R2 ( I , ISEG , IEDGE) ) 

2 /SIN(THETA(ISEG, IEDGE) ) 

210 TANGENT ( I ) = THETA( ISEG , IEDGE ) * 

1 (COS(PHI2)*R2(I, ISEG, IEDGE) - COS ( PHI1 ) *R1 ( I , ISEG , IEDGE ) ) 

2 /SIN(THETA( ISEG, IEDGE) ) 

GO TO 1000 


EDGE OF REVOLUTION MAP = 4 


PROJECTION ALONG AXIS AND RADIUS 

300 RP = ( RPA2 ( ISEG , IEDGE )-RPAl ( ISEG , IEDGE ) ) * RAT IO + RPA1 ( ISEG, IEDGE) 
RM = ( RM2( ISEG, IEDGE)- RM1 ( ISEG, IEDGE) ) * RATIO + RM1 { ISEG, IEDGE) 

— ANGLE 

GAMMA = THETA ( ISEG , IEDGE ) * RATIO 

— CALCULATE THE POSITION AND TANGENT 

DO 410 1=1,3 
C 
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UR = COS (GAMMA) *UI(I,ISEG,IEDGE) 

& + SIN( GAMMA) *UJ ( I , ISEG , I EDGE ) 

C 

POINT(I) = COEFE(I, ISEG, IEDGE) + RP*UK( I , ISEG, IEDGE) + RM*UR 
C 

410 TANGENT ( I ) = COS < GAMMA) *UJ ( I , ISEG , IEDGE ) 

& - SIN(GAMMA)*UI( I, ISEG, IEDGE) 

C 

GO TO 1000 


HGM HOLE MAP = 3 

400 IDUCT = IZONE - 3 
IF ( I DUCT . EQ . 1 ) THEN 
UPPER DUCT 


IF( IEDGE. EQ. 5) THEN 

ICLOCK = 0 

ANGLE = (1.0 - RATIO) *130.0 + RATIO*225.0 
END IF 


IF( IEDGE. EQ.6) THEN 


ICLOCK 

ANGLE 

END IF 

IF( IEDGE.EQ. 1 ) THEN 

ICLOCK 

ANGLE 

END IF 

IF( IEDGE.EQ. 9) THEN 

ICLOCK 

ANGLE 


1 

( 1.0 


- RATIO) *50.0 - RATIO* 50 . 0 


1 

( 1.0 


- RATIO) *130.0 + RATIO* 50 . 0 


0 

( 1.0 


- RATIO) *225.0 + RATIO*310.0 


END IF 


ELSE 


— LOWER DUCT 


IF( IEDGE.EQ. 1) 

THEN 


ICLOCK 

ANGLE 


END IF 

IF ( IEDGE.EQ. 5) 

THEN 


1 

(1.0 - RATIO) *130.0 ,+ RATIO* 50 . 0 
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C 

C 

C 

C 

C 

C 

C 


ICLOCK = 0 

ANGLE = (1.0 - RATIO) *130.0 + RATIO*180.0 
END IF 

IF ( IEDGE . BQ . 6 ) THEN 

ICLOCK = 1 

ANGLE = (1.0 - RATIO) *50.0 
END IF 

END IF 

CALL HOLES( ANGLE, POINT, TANGENT, ICLOCK , IDUCT) 

GO TO 1000 

HGM FAIRING MAP = 

500 CALL FAIRING( PT1 ( 1 , ISEG, IEDGE) , PT2 ( 1 , ISEG, IEDGE) , POINT, TANGENT, C) 

GO TO 1000 


9 


1000 CONTINUE 

RETURN 

END. 

SUBROUTINE SURFACE ( INIT , IS IDE ) 


THIS ROUTINE DETERMINES THE COORDINATES, FLOW VECTOR AND NORMAL OF A 
POINT ON A SURFACE. 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 


COMMON /COEFF/ 
COMMON /COORD/ 
COMMON /COUNTER/ 
COMMON /HEADER/ 
COMMON /HGM/ 
COMMON /HOTGAS/ 

COMMON /INITA/ 
COMMON /INITB/ 
COMMON /INITC/ 
COMMON /INPUT A/ 
COMMON /IWRITE/ 
COMMON /MARCHS/ 
COMMON /MAPED/ 
COMMON /MAPING/ 
COMMON /MAXIMUM/ 
COMMON /PARTIAL/ 
COMMON /SPACING/ 
COMMON /UNITS/ 


COEFE( 8,10,12) , COEFS (8,6), NMBRSEG( 12 ) 

UA( 3, 6), UE (3, 6), UN (3, 6) 

NODESAV , MATES AV , NODETOT , NBNODES , NPLANE 

ITITLE( 20) , LINE 

IHOLE 

IHGM, ETA1 1 , ETA1 2 , ETAl 3 , ETA3 1 , ETA32 , ETA33 , 

STR1 1 , STR1 2 , STR1 3 , STR3 1 , STR32 , STR33 , STR34 
IDIM, I FORMAT , MATRIX , MATING , MAPTEN , INCHES 
IBOX( 4,6), I RAY (4,3), NRAY ,IPT1(12) ,IPT2(12) 

PI , RADDEG 

EDGE( 6,12) , POINT{ 6 , 8 ) , SIDE( 6 , 6 ) , INPUTBC ( 6 ) 

IWRTC , IWRTI , IWRTN 
MARCH, INDEX( 3 ) 

KSEG ( 1 2 ) , UAXIS ( 3,6) ,IEDGEl( 6) ,IEDGE2(6) , GAMMA 
MAPSIDE ( 6 ) , MAPSEG( 10,12) 

ETAMAXf 10, 12) , SEGMAX ( 6 , 10 , 1 2 ) 

DEDN(3, 12) ,DSDN( 3, 2) ,SNORMAL( 3,6) 

ISTRTCH ( 3 ) , STRETCH ( 3 ) , ETAS ( 3 , 200 ) , ETA( 3 ) , DETA( 3 ) 
NU5 , NU6 ,NU19 , NU20 
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COMMON /ZONING/ ISECT, NSECT, IZONE , NMBRNDS ( 3 ) 

C 

DIMENSION EDG( 7 , 4 ) , EDGX{ 4 ) , EDGA( 4 ) , PT( 7 , 4 ) , PTX( 4 ) , PTA( 4 ) 
DIMENSION ET1 ( 4 ) ,ET3(4) ,SWE£P( 3) ,VDUM( 3) 

DIMENSION F(8),E1(3) ,E2( 3) ,R1 ( 3) ,R2( 3) 

DIMENSION ZERO(3),V£CTER(3) 

DIMENSION IETA1 ( 6 ) ,IETA2(6) ,SIGNS(6) 

DIMENSION UNITRAD( 3 ) 

C 

DATA IETA1 /l, 1,1, 1,2, 2/ 

DATA IETA2 /2 , 3 , 2 , 3 , 3 , 3/ 

DATA SIGNS /I . 0 , -1 . 0 , -1 . 0 , 1 . 0 , 1 . 0 ,-l . 0/ 

DATA ZERO /0 . ,0. ,0./ 

C 

NODNUM = NODESAV NODETOT + 1 
ETA COEFFICIENTS FOR BI-LINEAR INTERPOLATION 


F(l) = 1.0 - ETA( IETA1 ( IS IDE) ) 

F ( 2 ) = ETA( IETA1 ( ISIDE ) ) 

F(3) = 1.0 - ETA(IETA2( ISIDE)) 

P ( 4 ) = ETA( IETA2( ISIDE) ) 

F(5) = F ( 1 ) *F( 3 ) 

F(6) = F ( 1 ) *F ( 4 ) 

F( 7 ) = F( 2 ) *F ( 3 ) 

F(8) = F( 2) *F( 4 ) 

DETERMINE THE EDGES OF THE SURFACE 

LINE1 = IBOX(l, ISIDE) 

LINE2 = IBOX(2, ISIDE) 

LINE3 = IBOX(3, ISIDE) 

LINE4 = IBOX(4, ISIDE) 

DETERMINE THE CORNER POINTS OF THE SURFACE 

LPT 1 = IPTl(LINEl) 

LPT2 = IPT1 ( LINE3 ) 

LPT3 = IPT2 ( LINE3) 

LPT 4 = IPT2 ( LINE1 ) 


MAP = MAPSIDE( ISIDE) 

GO TO (100,100,300,400,500,600,700,800) MAP 
FLAT OR CYLINDRICAL SURFACE MAP = 1 OR 2 


CALCULATE POSITION 


C 


C 


100 DO 110 J=1 , 3 

SIDE(J, ISIDE) 

1 

2 

3 

DSDN ( J , 1 ) 

1 

2 

3 


= F(1)*EDGE(J, LINED 
+ F ( 2 ) ‘EDGE ( J , LINE3 ) 
+ F( 3 ) *EDGE( J , LINE2 ) 
+ F( 4 ) *EDGE ( J , LINE4 ) 

= EDGE ( J , LINE3 ) 

- EDGE(J, LINED 
+ F( 3 ) *DEDN ( J , LINE2 ) 
+ F( 4 ) *DEDN{ J , LINE4 ) 


- F ( 5 ) * POI NT ( J , LPT 1 ) 

- F( 6 ) * POINT ( J , LPT4 ) 

- F( 7 ) *POINT( J , LPT2 ) 

- F ( 8 ) *POINT( J , LPT3 ) 

+ F ( 3 ) * POI NT ( J , LPT 1 ) 
+ F( 4 ) *POINT( J , LPT4 ) 

- F( 3 ) *POINT( J , LPT2 ) 

- F ( 4 ) * PO I NT ( J , LPT3 ) 
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110 DSDN ( J / 2 ) = EDGE { J , LINE4 ) + F ( 1 ) * POINT* J , LPT1 ) 

1 - EDGE( J,LINE2) + F( 2 ) *P0INT( J , LPT2 ) 

2 + F( 1)*DEDN(J, LINED - F( 1 ) * P0INT( J , LPT4 ) 

3 + F<2)*DEDN(J,LINE3) - F ( 2 ) *P0INT( J , LPT3 ) 

C 

CALL CROSS {DSDN(1,1),DSDN( 1,2) ,SN0RMAL( 1 , ISIDE) ,40) 

C 

GO TO 1000 

SPECIAL SURFACE MAP = 3 

300 GO TO 1000 

EDGE OF REVOLUTION MAP = 4 

400 IF(INIT.EQ.l) THEN 
ETA DIRECTION OF REVOLUTION 

IONETWO = 1 

IF( COEFS ( 7 , ISIDE) . EQ . IETA2 (ISIDE ) ) IONETWO = 2 
EDGE 1 

IEDGEK ISIDE) = IBOX( IONETWO, ISIDE) 

EDGE 2 

IEDGE2( ISIDE) = IBOX( IONETWO + 2, ISIDE) 

UA: NORMALIZED VECTOR ALONG AXIS OF REVOLUTION 

CALL VADD( 1 . 0 , COEFS ( 4 , ISI DE ) , 1 . 0 , ZERO, VECTER, UA( 1 , ISIDE) ) 

UN: NORMALIZED VECTOR FROM AXIS TOWARD SURFACE 

CALL VADD( 1.0, EDGE( 1, IEDGEK ISIDE) ) ,- 1 . 0 , COEFS ( 1 , ISIDE ), E 1 , VECTER ) 
CALL CROSS (UA(1, ISIDE) , El , UN ( 1 , ISIDE ) , 4 1 ) 

UE: NORMALIZED VECTOR FROM AXIS TO FIRST EDGE 

CALL CROSS ( UN( 1, IS IDE) , UA( 1 , ISIDE ) , UE ( 1 , IS IDE ) , 4 2 ) 

END IF 

Rl: VECTOR FROM AXIS TO FIRST EDGE 

CALL VADD( 1.0, EDGE ( 1, IEDGEK ISIDE) ) , -1 . 0 , COEFS ( 1 , IS IDE ) , El, VECTER) 
CALL VDOT( El ,UA( 1 , ISIDE) , EPA1 ) 

CALL VADD(1.0, El, -EPA1,UA(1, ISIDE) ,R1, VECTER) 

CALL VMAG ( R 1 , RM 1 ) 

— DERIVATIVE 

CALL VDOT(DEDN( 1, IEDGEK ISIDE) ) ,UA( 1, ISIDE) ,DA1) 
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C 

CALL VDOT( DEDN{ 1 , I EDGE 1 ( ISIDE) ) , VECTER, DR1 ) 

C 

CALL VMAG( DEDN ( 1 , IEDGE1 ( ISIDE) ) , DMl ) 

— R2 : VECTOR FROM AXIS TO SECOND EDGE 

CALL VADD( 1 . 0 , EDGE ( 1 , IEDGE2 ( IS IDE ) ) ,-l .0,COEFS( 1 , ISIDE) ,E2, VECTER) 
CALL VDOT (E2,UA(1, ISIDE), EPA2 ) 

CALL VADD(1.0,E2,-EPA2,UA(1, ISIDE) ,R2, VECTER) 

CALL VMAG(R2,RM2) 

DERIVATIVE 

CALL VDOT ( DEDN ( 1, IEDGE2( ISIDE) ) , UA( 1 , IS I DE ) ,DA2) 

CALL VDOT ( DEDN ( 1, IEDGE2( ISIDE) ) , VECTER, DR2) 

CALL VMAG( DEDN( 1 , IEDGE2 ( ISIDE ) ) ,DM2) 

SWEEP ANGLE 

CALL VDOT ( R1 , R2 , R1DOTR2 ) 

ARG = R1DOTR2/C RM1*RM2 ) 

IF(ABS(ARG) .GT.1.0) ARG * ARG/ABS ( ARG ) 

THETA = ACOS(ARG) 

—CALCULATE RADIUS 

N = COEFS(7, ISIDE) 

— POSITION 

EP = ( EPA2 - EPAl ) *ETA( N) + EPA1 
RM = ( RM2 - RM1 ) * ETA ( N ) + RM1 

DERIVATIVE 

DA = ( DA2 - DAI ) *ETA( N) + DAI 
DR = ( DR2 - DR1 ) *ETA( N) + DR1 
DM = (DM2 - DMl ) *ETA( N) + DMl 


C 

C 

C 


ANGLE 

GAMMA = THETA*ETA( N) 

— CALCULATE THE POSITION AND SURFACE NORMAL 

DO 410 1=1,3 

UR = COS ( GAMMA) *UE ( I , IS I DE) + SIN( GAMMA) *UN( I , ISIDE) 
SIDEtl, ISIDE) = COEFS(I, ISIDE) + EP*UA( I , ISIDE) + RM*UR 
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DSDN ( I / 1 ) = { DR*UR + DA*UA< I , ISI DE ) )/DM 

C 

410 DSDN (1,2) = COS ( GAMMA ) * UN ( I , I S I DE ) - SIN( GAMMA) *UE( I , ISIDE ) 

— ORIENT SURFACE NORMAL DEPENDING ON DIRECTION OF EDGE REVOLUTION 
IF(COEFS(7, ISIDE) . EQ . IETA2 ( IS IDE ) ! THEN 
CALL CROSS ( DSDN ( 1, 1) , DSDN ( 1,2) ,SNORMAL( 1, IS IDE) ,43) 

ELSE 

CALL CROSS ( DSDN (1,2) , DSDN (1,1), SNORMAL( 1 , ISIDE) , 44) 

END IF 

GO TO 1000 


BOWL SURFACE WITH HOLES MAP = 5 


500 ETA1 = F( 2) 

ET1(1) = (ETA11 - 1.0) / REAL ( NMBRNDS ( 1 ) - 1) 

ET1 ( 2 ) = ( ETA 12 - 1.0) / REAL ( NMBRNDS ( 1 ) - 1) 

ET 1 ( 3 ) = ( ETA13 - 1.0) / REAL( NMBRNDS ( 1 ) - 1) 

ET 1 ( 4 ) = 1.0 

IET1 = 1 

IF(ETA1.GE.ET1(1) ) IET1 = 2 
IF ( ETAl . GE . ET1 ( 2 ) ) IET1 = 3 
IF( ETA1 .GE . ET1 ( 3 ) ) IET1 = 4 

ETA3 = F ( 4 ) 

ET3 ( 1 ) = ( ETA31 - 1.0) / REAL( NMBRNDS ( 3 ) - 1) 

ET3 ( 2 ) = ( ETA32 - 1.0) / REAL( NMBRNDS ( 3 ) - 1) 

ET3( 3 ) = ( ETA33 - 1.0) / REAL( NMBRNDS( 3) - 1) 

ET3 ( 4 ) = 1.0 

IET3 = 1 

I F ( ETA3 . GE . ET3 ( 1 ) ) IET3 = 2 
IF( ETA3.GE.ET3( 2) ) IET3 = 3 
IF( ETA3 .GE . ET3 ( 3 ) ) IET3 = 4 

IHOLE = 0 
ICLOCK = 0 

GO TO (510,520,540,560) IET1 


IET1 = 1 


510 CALL CONVERT ( EDGE( 1 , 3 ) , EDGX ( 1 ) , EDGA ( 1 ) ) 

CALL CONVERT { EDGE ( 1, 11) , EDGX ( 3 ) , EDGA( 3 ) ) 

RANGE = EDGX(l) 

ANGLE = (1.0 - ETA3) *EDGA( 1 ) + ETA3*EDGA(3) 
GO TO 580 


IET1 = 2 
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C 

520 RATI01 = { ETA1 - ET1 ( 1 ) ) / (ET1{2) - ET1 ( 1 ) ) 

C 

XI = RATIO 1*STR1 1/2.0 

ETAX1 = (EXP(Xl) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 

C 

X2 = STR11/2.0 

ETAX2 = (EXP(X2) - EXP( -X2 ) )/( EXP( X2 ) + EXP(-X2)) 

C 

EPS 1 = ETAX1/ETAX2 
C 

GO TO (522,524,526,528) IET3 

— IET3 = 1 

522 RATI03 = ETA3/ET3 ( I ) 

XI = RATIO3*STR31/2.0 

ETAX1 - (EXP(Xl) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 

X2 * STR31/2 . 0 

ETAX2 = (EXP(X2) - EXP(-X2) )/( EXP(X2) + EXP(-X2)) V 
EPS 3 = ETAX1/ETAX2 

CALL CONVERT ( EDGE ( 1,3) ,EDGX(1) , EDGA( 1 ) ) 

CALL CONVERT ( EDGE (1,8) ,EDGX(4) , EDGA{ 4 ) ) 

EDGX ( 4 ) * SEGMAX (1,1,3) - 6.07 

I DUCT = 1 
PT ( 7 , 3 ) = 130.0 

CALL HOLES ( PT( 7 , 3 ) , PT( 1 , 3 ) , VDUM , ICLOCK , IDUCT ) 

CALL CONVERT (PT( 1,3) ,PTX(3) ,PTA(3)) 

EDGX( 3 ) = (1.0 - EPS1 ) *EDGX ( 4 ) + EPS1*PTX(3) 

RANGE = (1.0 - EPS 3 ) *EDGX( 1 ) + EPS3*EDGX(3) 

EDGA( 2 ) = EPS3*PTA( 3 ) 

ANGLE = (1.0 - EPS1 ) *EDGA( 4 ) + EPS1*EDGA(2) 

GO TO 580 

— IET3 = 2 


524 

C 

c 

c 

c 


RATI03 

XI 

ETAMID 

X2 

ETAMX 

EPS3 


( ETA3 - ET3 ( 1 ) )/ ( ET3 ( 2 ) - ET3( 1 ) ) 

( STR32/2 . 0 )/2 . 0 

(EXP(Xl) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 
( RATI03 - 0.5) *STR32/2 . 0 

(EXP(X2) - EXP(-X2 ) )/{ EXP( X2 ) + EXP(-X2)) 
(ETAMX + ETAM I D ) / ( 2 . 0 * ETAM ID) 
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C 

C 

C 

C 

C 

C 

C 


CALL CONVERT { EDGE( 1,8), EDGX( 4 ) , EDGA( 4 ) ) 

EDGX( 4) = SEGMAX( 1,1,3) - 6.07 
IDUCT = 1 

EDG( 7,2) = (1.0 - EPS3 ) *130.0 + EPS3*225.0 
CALL HOLES(EDG(7,2) , EDG( 1 , 2 ) , VDUM , ICLOCK , IDUCT ) 

CALL CONVERT ( EDG ( 1 , 2 ) , EDGX( 2 ) , EDGA( 2 ) ) 

RANGE = (1.0 - EPS 1 ) *EDGX( 4 ) + EPS1*EDGX(2) 

ANGLE = (1.0 - EPS 1 ) *EDGA( 4 ) + EPS1*EDGA(2) 

GO TO 580 

• — IET3 = 3 

526 RATI03 = ( ETA3 - ET3 ( 2 ) ) / ( ET3 ( 3 ) - ET3(2)) 

XI = ( STR33/2 . 0 )/2 . 0 

ETAMID = (EXP(Xl) - EXP( -Xl ) ) /( EXP( XI ) + EXP(-Xl)) 

X2 = ( RATI03 - 0.5)*STR33/2.0 
ETAMX = (EXP(X2) - EXP ( -X2 ) ) / ( EXP( X2 ) + EXP(-X2)) 

EPS3 = (ETAMX + ETAMID)/(2.0*ETAMID) 

CALL CONVERT (EDGE (1,8) ,EDGX(4) ,EDGA{4) ) 

EDGX ( 4 ) = SEGMAX (1,1,3) - 6.07 

IDUCT = 1 
PT ( 7 , 2 ) = 225.0 

CALL HOLES ( PT( 7, 2) ,PT( 1,2) ,VDUM, ICLOCK, IDUCT) 

CALL CONVERT ( PT ( 1 , 2 ) , PTX( 2 ) , PTA( 2 ) ) 

IDUCT = 2 
PT ( 7 , 3 ) = 130.0 

CALL HOLES ( PT( 7, 3 ) , PT( 1,3), VDUM , ICLOCK, IDUCT) 

Call convert ( pt( l , 3 ) , ptx{ 3 ) , pta( 3 ) ) 

EDGX( 2 ) = (1.0 - EPS3 ) *PTX( 2 ) + EPS3*PTX(3) 

RANGE = (1.0 - EPS1 ) *EDGX( 4 ) + EPS1*EDGX(2) 

E DGA ( 2 ) = (1.0 - EPS3) *PTA{ 2) + EPS3*PTA(3) 

ANGLE = (1.0 - EPS1 ) *EDGA( 4 ) + EPS1*EDGA(2) 

GO TO 580 

— IET3 = 4 

528 RATI03 = 1.0 - ( ETA3 - ET3 ( 3 ) ) /( ET3 ( 4 ) - ET3(3)) 
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C 

C. 

C 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


XI = RATI03*STR34/2 • 0 

ETAX1 = (EXP(Xl) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 

X2 = STR34/2 . 0 

ETAX2 = (EXP(X2) - EXP(-X2 ) )/ ( EXP ( X2 ) + EXP(-X2)) 
EPS3 » 1.0 - ETAX1/ETAX2 
CALL CONVERT ( EDGE ( 1,8) ,EDGX(4) , EDGA( 4 ) ) 

EDGX( 4 ) = SEGMAX{ 1,1,3) - 6.07 
IDUCT = 2 

E DG ( 7 / 2 ) = (1.0 - EPS3 ) * 130 . 0 + EPS3*180.0 
CALL HOLES ( EDG( 7,2) , EDG( 1 , 2 ), VDUM , ICLOCK , IDUCT) 

CALL CONVERT (EDG( 1,2) ,EDGX(2) , EDGA( 2 ) ) 

RANGE = (1.0 - EPS1 ) *EDGX( 4 ) + EPS1*EDGX(2) 

ANGLE = (1.0 - EPS1 ) *EDGA( 4 ) + EPS1*EDGA(2) 

GO TO 580 
IET1 = 3 

540 RATIOl = ( ETA1 - ET1 ( 2 ) )/( ET1 ( 3 ) - ET1(2)) 

XI = ( STR1 2/2 . 0 )/2 . 0 

ETAMID = ( EXP ( XI ) - EXP( -XI ) ) / ( EXP( XI ) + EXP(-Xl)) 

X2 = (RATIOl - 0.5) ‘STR12/2 . 0 
ETAMX » (EXP(X2) - EXP( -X2 ) ) / ( EXP( X2 ) + EXP(-X2)) 

EPS1 = (ETAMX + ETAMID) /( 2.0* ETAMID) 

GO TO (542,544,546,548) IET3 

— IET3 = 1 

542 RATI03 = ETA3/ET3(1) 

XI = RATI03*STR3 1/2.0 

ETAX1 = (EXP(Xl) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 
X2 = STR31/2.0 

ETAX2 = (EXP(X2) - EXP { -X2 ) ) /( EXP( X2 ) + EXP(-X2)) 
EPS3 = ETAX1/ETAX2 

CALL CONVERT ( EDGE ( 1 , 3 ) , EDGX ( 1 ) , EDGA ( 1 ) ) 

IDUCT = 1 

EDG( 7,3) = (1.0 - EPS1 )* 130 .0 + EPS1*50.0 
CALL HOLES ( EDG( 7 , 3 ) , EDG( 1,3) , VDUM , ICLOCK , IDUCT) 

CALL CONVERT (EDG( 1,3) ,EDGX( 3) ,EDGA( 3) ) 


c 
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C 

C 

c 


RANGE = 

(1.0 

ANGLE = 

(1.0 

GO TO 580 



EPS3)*EDGX( 1) + EPS3*EDGX( 3 ) 
EPS3 ) *EDGA( 1 ) + EPS3*EDGA( 3 ) 


C 

C IET3 

C 

544 

C 


C 


C 

c 


RATIO 3 = 

XL = 
ETAMID = 

X2 = 
ETAMX . = 

EPS3 = 

IHOLE = 
I DUCT = 


C 


EDG( 7 / 1 ) = 
EDG( 7,2) = 
EDG( 7,3) = 
EDG( 7,4) = 

C 

PT ( 7 , 1 ) = 
PT (7,2) = 
PT ( 7 , 3 ) = 
PT (7,4) = 
C 


(ETA3 - ET3 ( 1 ) ) /( ET3 ( 2 ) - ET3 ( 1 ) ) 


( STR32/2 . 0 ) /2 . 0 

(EXP(Xl) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 


( RATI03 - 0.5) *STR32/2 .0 

(EXP(X2) - EXP ( -X2 ) )/( EXP( X2 ) + EXP(-X2)) 

(ETAMX + ETAMID)/ ( 2 .0* ETAMID) 

1 

1 


( 1.0 

( 1.0 

( 1.0 

( 1.0 


EPS1 ) * 130 . 0 + 
EPS3 ) *50.0 + 
EPS1 ) *225 . 0 + 
EPS3)*130.0 + 


EPS1*50.0 
EPS3*(-50.0) 
EPS 1*310.0 
EPS3* 225 . 0 


130.0 
50.0 

310.0 

225.0 


DO 545 N=1 , 4 
C 

CALL HOLES (EDG( 7, N) , EDG( 1 , N) , VDUM , ICLOCK , IDUCT ) 
C 


CALL CONVERT{ EDG( 1 , N) ,EDGX(N) , EDGA( N ) ) 

C 

CALL HOLES ( PT( 7 , N) ,PT(1,N) , VDUM , ICLOCK , IDUCT) 
C 


C 


c 


c 


545 CALL CONVERT (PT(1,N) ,PTX(N) ,PTA(N)) 


RANGE = (1.0 

1 

2 

3 + (1.0 


- EPS3 ) *EDGX( 1 ) 
+ EPS1*EDGX( 2 ) 
+ EPS3*EDGX( 3 ) 

- EPS1 ) *EDGX( 4 ) 


(1.0 - EPS1 ) * ( 1 . 0 - EPS 3 ) * PTX( 1 ) 

- (1.0 - EPS 3 ) * EPS 1 * PTX ( 2 ) 

- EPS 1 * EPS 3* PTX ( 3 ) 

- (1.0 - EPS1 ) * EPS 3* PTX ( 4 ) 


ANGLE = (1.0 

1 

2 

3 + (1.0 


- EPS3 ) *EDGA( 1 ) 
+ EPS1* EDGA( 2 ) 
+ EPS3*EDGA( 3 ) 

- EPS1 ) *EDGA( 4 ) 


(1.0 - EPS1 ) *(1.0 - EPS 3 ) *PTA( 1 ) 

- (1.0 - EPS 3 ) *EPS1* PTA( 2 ) 

- EPS1*EPS3*PTA( 3 ) 

- (1.0 - EPS1 ) *EPS 3* PTA( 4 ) 


GO TO 580 


C 

C IET3 = 3 

C 

546 RATIO 3 = ( ETA3 - ET3 ( 2 ) ) / ( ET3 ( 3 ) - ET3(2)) 

C 

XI = ( STR33/2 . 0 )/2 .0 

ETAMID = (EXP(Xl) - EXP( -XI ) ) /{ EXP( XI ) + EXP(-Xl)) 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


X2 = ( RAT 10 3 - 0.5) *STR33/2 . 0 
ETAMX = (EXP(X2) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 

EPS3 = (ETAMX + ETAMID) / ( 2 . 0*ETAMID) 

I DUCT = 1 

E DG ( 7 , 1 ) = (1.0 - EPS1 ) * 22 5 . 0 + EPS1*310.0 
.CALL HOLES ( EDG( 7 , 1 ) , EDG( 1 , 1) , VDUM , I CLOCK , I DUCT ) 

CALL CONVERT ( EDG( 1,1), EDGX( 1 ) , EDGA( 1 ) ) 

1DUCT = 2 

EDG( 7,3) = (1.0 - EPS 1 ) * 130 . 0 + EPS1*50.0 
CALL HOLES ( EDG ( 7 , 3 ) , EDG ( 1,3), VDUM , ICLOCK , IDUCT ) 

CALL CONVERT (EDG (1,3) ,EDGX(3) , EDGA( 3)) 

RANGE = (1.0 - EPS3 ) * EDGX( 1 ) + EPS3*EDGX(3) 

ANGLE = (1.0 - EPS 3 ) *EDGA( 1 ) + EPS3*EDGA(3) 

GO TO 580 

— IET3 = 4 


C 


C 


C 

C 


C 


C 


C 

C 

C 

C 

C 

C 


548 RATI03 = 1.0 - ( ETA3 - ET3 ( 3 ) )/( ET3 ( 4 ) - ET3 ( 3 ) ) 

XI = RATI03 *STR34/2 . 0 

ETAXl = (EXP(Xl) - EXP ( -XI ))/(EXP(Xl) + EXP(-Xl)) 
X2 = STR34/2 . 0 

ETAX2 = (EXP(X2) - EXP( -X2 ) )/( EXP( X2 ) + EXP(-X2)) 

EPS3 = 1.0 - ETAX1/ETAX2 

IHOLE = 1 
IDUCT = 2 

E DG ( 7 , 1 ) = (1.0 - EPS1 ) * 130 . 0 + EPS1*50.0 

EDG( 7,2) = (1.0 - EPS3 ) *50 . 0 

EDG( 7,4) = (1.0 - EPS3) * 130 .0 + EPS3*180.0 

PT (7,1) = 130.0 
PT ( 7 , 2 ) = 50.0 

PT ( 7 , 3 ) = 0.0 

PT( 7 , 4 ) = 180.0 

DO 549 N=1 , 4 

CALL HOLES (EDG( 7, N) ,EDG( 1,N) ,VDUM, ICLOCK, IDUCT) 

CA^.L CONVERT ( EDG ( 1,N) , EDGX( N) , EDGA( N ) ) 

CALL HOLES ( PT( 7 , N) , PT( 1 , N) , VDUM , ICLOCK , IDUCT ) 

549 CALL CONVERT ( PT( 1 , N) , PTX( N) , PTA( N) ) 


C-35 

LOCKHtED-HUNTSVILLE ENGINEERING CENTER 



non on o n non 


LMSC-HEC TR F225755 


C 

C 


C 


EDGX( 3 ) 
EDGA( 3) 


( 1.0 

( 1.0 


EPS 1 ) *PTX( 4 ) + EPS1*PTX(3) 
EPS1 ) *PTA( 4 ) + EPS1*PTA( 3) 


RANGE = (1.0 

1 

2 

3 + ( 1.0 


- EPS3 ) * EDGX( 1 ) 
+ EPS1*EDGX( 2) 
+ EPS3*EDGX(3) 

- EPS1 ) *EDGX( 4 ) 


(1.0 - EPS1 ) * ( 1 . 0 - EPS3 ) * PTX( 1 ) 

- (1.0 - EPS3 ) *EPS 1*PTX( 2 ) 

- EPS1*EPS3*PTX( 3 ) 

- (1.0 - EPS1 ) *EPS3*PTX( 4 ) 


ANGLE = (1.0 

1 

2 

3 + (1.0 


- EPS3 ) *EDGA( 1 ) 
+ EPS 1*EDGA( 2 ) 
+ EPS3*EDGA( 3 ) 

- EPS1 ) *EDGA( 4 ) 


(1.0 - EPS 1 ) * ( 1 . 0 - EPS 3 ) *PTA( 1 ) 

- (1.0 - EPS3 ) *EPS1*PTA( 2 ) 

- EPS1*EPS3*PTA( 3) 

- (1.0 - EPS1 ) *EPS3*PTA( 4 ) 


GO TO 580 


IET1 = 4 


560 RATIOl 
XI 

ETAMID 

X2 

ETAMX 

EPS 1 

GO TO ( 

IET3 = 1~ 


= ( ETA1 - ET1 ( 3 ) )/ ( ET1 ( 4 ) - ET1(3)) 

= ( STR13/2 . 0 )/2 . 0 

= ( EXP ( XI ) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 
= (RATIOl - 0.5) *STR13/2 . 0 

= ( EXP(X2 ) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 
= (ETAMX + ETAMID)/( 2 • 0* ETAMID) 
562,564,566,568) IET3 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


562 


RATI03 = ETA3/ET3 ( 1 ) 

XI = RAT 103 *STR 3 1/2.0 
ETAX1 = (EXP(Xl) - EXP ( -XI ) ) /( EXP ( XI ) + 

X2 = STR31/2 . 0 

ETAX2 = (EXP(X2) - EXP(-X2) )/( EXP(X2 ) + 
EPS 3 = ETAX1/ETAX2 

CALL CONVERT ( EDGE( 1,3), EDGX( 1 ) , EDGA( 1 ) ) 

CALL CONVERT ( EDGE { 1 , 7 ) , EDGX( 2 ) , EDGA( 2 ) ) 

I DUCT = 1 
PT ( 7 , 4 ) = 50.0 

CALL HOLES ( PT( 7 , 4 ) , PT( 1 , 4 ) , VDUM, ICLOCK , IDUCT) 
CALL CONVERT ( PT( 1,4), PTX( 4 ) , PTA( 4 ) ) 

EDGX( 3 ) = (1.0 - EPS1 ) * PTX( 4 ) + EPS1*EDGX( 
RANGE = (1.0 - EPS3 ) *EDGX( 1 ) + EPS3*EDGX( 
EDGA( 4 ) = EPS3*PTA( 4 ) 

ANGLE = (1.0 - EPS1 ) *EDGA( 4 ) + EPS1*EDGA( 


EXP(-Xl)) 
EXP(-X2) ) 


2 ) 

3) 

2 ) 
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C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


GO TO 580 

— IET3 » 2 

564 RATI03 = ( ETA3 - ET3( 1 ) )/( ET3( 2) - ET3 ( 1 ) ) 

XI = ( STR32/2 . 0 )/2 . 0 

ETAMID = ( EXP ( XI ) - EXP ( -XI ) ) / ( EXP ( XI ) + EXP(-Xl)) 

X2 = ( RATI03 - 0.5) ‘STR32/2 . 0 
ETAMX = (EXP(X2) - EXP(-X2) )/(EXP( X2) + EXP(-X2)) 

EPS3 = (ETAMX + ETAMID) / (2 . 0 ‘ETAMID ) 

CALL CONVERT ( EDGE( 1 / 7 ) , EDGX( 2 ) , EDGA{ 2) ) 

IDUCT = 1 

EDG( 7,4) = (1.0 - EPS3 ) * 50 . 0 + EPS3*(-50.0) 
CALL HOLES (EDG(7,4) ,EDG(1,4) , VDUM , IC LOCK , IDUCT) 

CALL CONVERT ( EDG( 1,4), EDGX( 4 ) , EDGA( 4 ) ) 

RANGE = (1.0 - EPS1 ) *EDGX( 4 ) + EPS1*EDGX(2) 

ANGLE = (1.0 - EPS1 ) *EDGA( 4 ) + EPS1*EDGA(2) 

GO TO 580 

— IET3 = 3 

566 RATI03 = ( ETA3 - ET3( 2) )/(ET3( 3) - ET3(2)) 

XI = ( STR33/2 . 0 )/2 . 0 

ETAMID = ( EXP ( XI ) - EXP( -XI ) ) / ( EXP( XI ) + EXP(-Xl)) 

X2 = ( RAT IO 3 - 0.5) ‘STR33/2 . 0 
ETAMX =* (EXP(X2) - EXP( -X2 ) ) / ( EXP( X2 ) + EXP(-X2)) 

EPS3 = (ETAMX + ETAM I D ) / ( 2 . 0 * ETAM ID) 

CALL CONVERT ( EDGE( 1,7), EDGX( 2 ) , EDGA( 2 ) ) 

IDUCT = 1 
PT ( 7 , 1 ) = 310.0 

CALL HOLES ( PT ( 7 , 1 ) , PT ( 1 , 1 ) , VDUM , I CLOC K , I DUCT ) 

CALL CON VE RT ( PT ( 1 , 1 ) ,PTX(1) ,PTA(1) ) 

IDUCT = 2 
PT ( 7 , 4 ) = 50.0 

CALL HOLES ( PT( 7 , 4 ) , PT( 1 , 4 ) , VDUM , ICLOCK , IDUCT) 

CALL CONVERT ( PT( 1,4), PTX( 4 ) , PTA( 4 ) ) 

EDGX( 4 ) = (1.0 - EPS 3 ) ‘ PTX( 1 ) + EPS3»PTX(4) 


C-37 

LOCKflIED-HUNTSVILLE ENGINEERING CENTER 



non o o onnnononoonnnn non n noon 


LMSC-HEC TR F225755 


RANGE = (1.0 - EPS1 ) *EDGX( 4 ) + EPS1*EDGX(2) 
EDGA( 4 ) = (1.0 - EPS3)*PTA( 1) + EPS3*PTA(4) 

ANGLE = (1.0 - EPS1) *EDGA( 4 ) + EPS1*EDGA(2) 

GO TO 580 

IET3 = 4 

568 RATI03 = 1.0 - ( ETA 3 - ET3( 3 ) )/( ET3 ( 4 ) - ET3{3)) 

XI = RATI03*STR34/2 . 0 

ETAX1 = ( EXP ( XI ) - EXP ( -XI ) )/( EXP (XI ) + EXP(-Xl)) 
X2 = STR34/2 . 0 

ETAX2 = (EXP(X2) - EXP( -X2 ) )/( EXP( X2 ) + EXP(-X2)) 
EPS3 = 1.0 - ETAX1/ETAX2 
CALL CONVERT ( EDGE (1,7), EDGX ( 2 ) , EDGA ( 2 ) ) 

IDUCT = 2 

EDG( 7,4) = (1.0 - EPS3)*50 . 0 
CALL HOLES (EDG (7,4) ,EDG( 1,4) , VDUM , ICLOCK , IDUCT) 

CALL CONVERT ( EDG( 1 , 4 ) , EDGX( 4 ) , EDGA( 4 ) ) 

RANGE = (1.0 - EPS1 ) *EDGX( 4 ) + EPS1*EDGX(2) 
ANGLE = (1.0 - EPS1 ) *EDGA( 4 ) + EPS1*EDGA(2) 

GO TO 580 
580 CONTINUE 

RADIUS AND TANGENT 


XX = RANGE + 6.07 

— FIRST SEGMENT 

I F ( XX . LE . S EGMAX (1,1,3)) THEN 

RAD = 7.434 

DSDN (1,1) = 1.0 
DSDN( 2,1) = 0.0 
DSDN( 3 , 1 ) = 0.0 

GO TO 590 

END IF 

— SECOND SEGMENT 

IF(XX.LT.SEGMAX( 1,2, 3) ) THEN 
C 

RAD = SQRT(7. 890**2 - (3.375 - RANGE) **2) 
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C 


C 


C 

C 


— LAST SEGMENT 


POSITION 


TH = ACOS( RAD/7. 890) 
IF( XX. GT .9.445) TH = -TH 

DSDN ( 1 , 1 ) = COS(TH) 

DSDN(2,1) = S IN{ TH ) 

DSDN (3,1) = 0.0 

GO TO 590 

END IF 


RAD = 14.5636605 - 0 . 57735027*XX 
TH = -30 . 0*RADDEG 


DSDN (1,1) = COS(TH) 
DSDN ( 2 , 1 ) = SIN(TH) 
DSDN ( 3 , 1 ) = 0.0 


590 SWEEP(l) = 0.0 

SWEEP ( 2 ) = COS (ANGLE) 

SWEEP( 3 ) = SIN( ANGLE) 

UA( 1 , ISIDE ) * 1.0 
UA( 2 , IS IDE) = 0.0 
UA(3, ISIDE) = 0.0 

CALL VADD( XX, UA( 1 , ISIDE ), RAD, SWEEP ,S IDE ( 1 , ISIDE ) ,VDUM) 
— NORMAL — — 


TANGENT 1 


DA = DSDN (1,1) 

DN = DSDN( 2 , 1 ) 

CALL VADD ( DN , SWEEP , DA , UA ( 1 , I S I DE ) , VDUM , DSDN ( 1 , 1 ) ) 

■ — TANGENT 2 

DSDN (1,2) = 0.0 

DSDN (2,2) = -SIN( ANGLE) 

DSDN( 3, 2) = COS ( ANGLE ) 

— SURFACE NORMAL 

CALL CROSS ( DSDN ( 1,1), DSDN (1,2) ,SNORMAL( 1, ISIDE) ,45) 
GO TO 1030 


MATED SURFACE OF DUCT MAP = 6 


600 EPS1 = F ( 2 ) 
EPS3 = F ( 4 ) 
C 
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C 


IDUCT = IZONE - 3 
IF(IDUCT.EQ.l) THEN 


— UPPER DUCT- 




* " “ 


EDG( 7,1) 

3 

(1.0 - 

EPS1 ) * 130 . 0 


EPS1*50.0 

EDG ( 7,2) 

s 

(1.0 - 

EPS3 ) * 50 . 0 


EPS3* (-50.0) 

EDG( 7,3) 

= 

(1.0 - 

EPS1 ) * 225 . 0 

♦ 

EPS1*310.0 

EDG (7,4) 

= 

(1.0 - 

EPS3 ) * 1 30 . 0 

+ 

EPS 3*225.0 

PT ( 7 , 1 ) 

s 

130.0 




PT ( 7 , 2 ) 

= 

50.0 




PT ( 7 » 3 ) 

= 

310.0 




PT ( 7 , 4 ) 

= 

225.0 





ELSE 


— LOWER DUCT- 


EDG( 7,1) = 

(1.0 - 

EPS1 ) * 130 . 0 + EPS1*50 . 0 

EDG ( 7,2) = 

(1.0 - 

EPS3 ) * 50 . 0 

EDG( 7,4) = 

(1.0 - 

EPS3)*130 .0 + EPS 3*180.0 

PT ( 7 , 1 ) = 

130.0 


PT( 7 , 2 ) = 

50.0 


PT ( 7 , 3 ) = 

0.0 


PT( 7 , 4 ) = 

180.0 



END IF 

CALCULATE POINTS AND EDGES 

DO 610 N=l,4 

CALL HOLES ( EDG( 7 , N) ,EDG(1,N) ,VDUM , ICLOCK , IDUCT ) 

CALL CONVERT ( EDG( 1 * N) , EDGX( N) , EDGA( N) ) 

CALL HOLES ( PT ( 7 , N ) , PT ( 1 , N ) ,VDUM , ICLOCK , IDUCT ) 

610 CALL CONVERT ( PT( 1 , N) , PTX( N) , PTA( N) ) 

IF( IDUCT. EQ. 2) THEN 

EDGX( 3 ) = (1.0 - EPS 1 ) *PTX( 4 ) + EPS1*PTX(3) 

EDGA( 3 ) = (1.0 - EPS 1 ) *PTA( 4 ) + EPS1*PTA(3) 

END IF 

— DISTANCE ALONG AXIS AND ANGLE FROM VERTICAL 

RANGE = (1.0 - EPS3) *EDGX( 1 ) - (1.0 - EPS1)*(1.0 - EPS3)*PTX(1) 


1 + EPS1*EDGX( 2 ) 

2 + EPS3*EDGX( 3) 

3 + (1.0 - EPS1 ) *EDGX( 4 ) 


- (1.0 - EPS 3) *EPS1*PTX( 2) 

- EPS1*EPS3*PTX( 3) 

- (1.0 - EPS1 ) * EPS3*PTX( 4 ) 


ANGLE = (1.0 - EPS 3 ) *EDGA( 1 ) - (1.0 - EPS1)*(1.0 - EPS3)*PTA(1) 
1 + EPS1*EDGA( 2 ) - (1.0 - EPS 3) *EPS1*PTA( 2) 
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2 + EPS3*EDGA( 3 ) - EPS1*EPS3*PTA( 3 ) 

3 + (1.0 - EPS1 ) *EDGA( 4 ) - (1.0 - EPS1 ) *EPS3*PTA( 4 ) 

— POSITION 


RAD = SQRT(7. 890**2 - ( 3.375 - RANGE ) ** 2 ) 

SIDE ( 1 , ISIDE ) = RANGE + 6.07 
SIDE( 2, IS IDE) = RAD*COS( ANGLE) 

SIDE(3, ISIDE) = RAD*S IN (ANGLE) 

— NORMAL 


VDUM(l) = SIDE( 1, ISIDE) - 9.445 
VDUM ( 2 ) = SIDE(2, ISIDE) 

VDUM ( 3 ) = SIDE( 3, ISIDE) 

CALL VMAG ( VDUM , SMAG ) 

SNORMALt 1, ISIDE) = VDUM(1)/SMAG 
SNORMAL( 2 , ISIDE ) = VDUM(2)/SMAG 
SNORMAL( 3, ISIDE) = VDUM(3)/SMAG 

GO TO 1030 


OUTER SURFACE OF DUCT MAP = 7 


700 IF ( ISIDE.LT. 4) THEN 
—SURFACE 1 AND 3 

M = LINE2 
N = LINE4 


ELSE 

— SURFACE 5 AND 6 

M = LINE1 
N = LINE3 

END IF 

CALL FAIRING ( EDGE ( 1, M) , EDGE (1,N) , S IDE ( 1 , IS IDE ) , 
* VDUM, SNORMAU 1, ISIDE ) ) 


GO TO 1030 


SPECIAL SURFACE MAP = 8 


800 GO TO 1000 


DIRECT SURFACE NORMAL INTO FLOW DOMAIN 


DO 1010 1=1,3 

1010 SNORMAL( I, ISIDE) = SNORMAL( I , ISIDE) *SIGNS( ISIDE) 


ERROR MESSAGE IF SURFACE NORMAL = 0 


1030 CALL VMAG(SNORMAL( 1, ISIDE) ,SNMAG) 
C 
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C 

C 

C 

C 

C 

C 


IF(SNMAG.EQ.O.O) THEN 

LINE = LINE + 5 
IF (LINE. GE. 60) THEN 

WRITE(NU6, 1100) I TITLE, ISECT,NSECT, I20NE 
LINE = 6 

END IF 

WRITE(NU6, 1110) 


RETURN 
END IF 

INTERMEDIATE PRINT 

I F ( IWRTI.EQ.O .OR. MOD( NODNUM , IWRTI ) . NE . 0 ) RETURN 
LINE = LINE + 3 
IF( LINE. GE. 60) THEN 

WRITE(NU6,1100) ITITLE, ISECT,NSECT,IZONE 
LINE = 4 
END IF 

RETURN 
FORMAT STATEMENTS 

1100 FORMAT ( 1H1 , 10X, 20A4 , 13X , 8H SECTION, 12 , 3H OF,I3,9H FOR ZONE, 13) 

1110 FORMAT (/ 49H SURFACE NORMAL EQUALS ZERO IN SUBROUTINE SURFACE) 

END 

*********************************************************************** 
• * ♦ * [ij* *********************** * ************ ******************* ****** * 

*********************************************************************** 

SUBROUTINE CONVERT ( POINT, RANGE , ANGLE ) 

CALCULATES PROJECTION AND SWEEP ANGLE OF POINT ON BOWL SURFACE 
POINT = POSITION ON BOWL SURFACE 

RANGE = PROJECTION ONTO X AXIS FROM EDGE OF BOWL 
ANGLE = SWEEP ANGLE FROM VERTICAL 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION POINT( 3 ) , BOWL( 3 ) , AXIS ( 3 ) , VECTOR ( 3 ) , SWEEP ( 3 ) , VDUM( 3 ) 
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C 

DATA BOWL /6 . 07 , 0 . 0 , 0 . 0/ 

DATA AXIS /l. 0,0. 0,0.0/ 

DISTANCE ALONG AXIS FROM EDGE OF BOWL 

ONE =1.0 

CALL VADD ( ONE , POI NT , -ONE , BOWL , VECTOR, VDUM ) 

CALL VDOT( VECTOR, AXIS, RANGE) 

ANGLE FROM VERTICAL 

CALL VADD( ONE, VECTOR, -RANGE, AXIS, VDUM, SWEEP) 

ANGLE = ACOS(SWEEP( 2) ) 

RETURN 

END 

*********************************************************************** 

* * * * * OUTPUT* *********************************************************** 
*********************************************************************** 

SUBROUTINE BLKOUT ( NUNIT , NODSTOR) 


WRITES THE FORMATTED BLOCKED GEOMETRY FILE (NUNIT) 


DOUBLE PRECISION NODE 
REAL NODE 

COMMON /INITA/ IDIM , MAPTEN, INCHES 

COMMON / IOCOUNT/ IREWIND( 40 ) , NREAD( 40 ) , NWRITE ( 40 ) 

COMMON /OUT/ NODE ( 5, 10000) 

COMMON /ZONING/ ISECT, NSECT, I ZONE , NMBRNDS ( 3 ) 

COMMON /MARCHS/ MARCH , INDEX( 3 ) 

IPLANE = NWRITE( NUNIT) + 1 

WRITE (NUNIT, 1000) NODSTOR, IPLANE ,( NMBRNDS ( I ) ,1=1,3) , MARCH 
WRITE (NUNIT, 10 10) ( NODE ( 1 , I ), 1=1 , NODSTOR) 

WRITE ( NUNIT ,1010) ( NODE ( 2 , I ), 1= 1 , NODSTOR) 

WRITE ( NUNIT ,1010) ( NODE ( 3 , I ), 1=1 , NODSTOR) 

RETURN 

FORMAT STATEMENTS 

1000 FORMAT ( 2415) 

1010 FORMAT ( 10E13 . 7 ) 

END 


SUBROUTINE REWRITE 
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THIS PROGRAM MAKES THE HGM GEOMETRY OUTPUT FILE 
COMPATABLE WITH PL0T3D GEOMETRY FILES 


COMMON /ZONING/ ISECT , NSECT, IZONE , NMBRNDS ( 3 ) 

COMMON /UNITS/ NU5 , NU6 , NU19 , NU20 , NU2 1 
COMMON /MARCHS/ MARCH , INDEX( 3 ) 

DIMENSION ID(10) ,JD(10) ,KD( 10) 

DIMENSION XBUF( 155000) , YBUF( 155000 ) , ZBUF ( 1 55000 ) 

REWIND ( NU20 ) 

12 = 0 
NGRID = 0 

READ( NU20 , 9000 , END= 200 ) NSTORE , IPLN, INOD2 , JNOD2 , KNOD2 , 
* MARCH 2 


9000 FORMAT ( 815) 

WRITE ( NU6 , 9000 ) NSTORE , IPLN , INOD2 , JNOD2 ,KNOD2 , MARCH2 

11 = 12 + 1 

12 = 12 + NSTORE 

READ(NU20, 9010, ERR= 400 ) < XBUF ( I) ,1=11,12) 

9010 FORMAT ( 10E13.7) 

READ( NU20 ,9010, ERR=400 ) ( YBUF ( I ) , 1 = 1 1 , 12 ) 

READ( NU 20, 9010, ERR=400 ) ( ZBUF( I ) , 1=11 , 12 ) 


100 READ( NU20, 9000 , END=200) NSTORE, IPLN, INOD, JNOD, KNOD, MARCH 
WRITE ( NU6 ,9000) NSTORE , IPLN , INOD , JNOD, KNOD , MARCH 


IF( INOD. NE.INOD2. OR. JNOD. NE.JNOD2 .OR. KNOD. NE .KNOD2 )THEN 
NGRID = NGRID + 1 
C 

IF ( MARCH 2 . EQ . 1 ) THEN 
ID(NGRID) = KNOD2 
JD( NGRID) = JNOD2 
KD( NGRID) = INOD2 
END IF 
C 

I F ( MARCH 2 . EQ . 2 ) THEN 
ID( NGRID) = INOD2 
JD( NGRID) = KNOD2 
KD( NGRID) = JNOD2 
END IF 
C 

IF ( MARCH 2 . EQ . 3 )THEN 
ID(NGRID) = JNOD2 
JD( NGRID) = INOD2 
KD{ NGRID) = KNOD2 
END IF 
C 
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WRITE ( NU6 (9120) NGRID, ID( NGRID) , JD( NGRID) ,KD(NGRID) ,12 
9120 FORMAT! ' NGRID, IDIM, JDIM, KDIM, I 2 ',515) 

C 

WRITE! 19,9010, ERR=450 ) 

1 (XBUF( I) ,1=1,12) , 

2 ( YBUF( I) , 1=1 , 12) , 

3 (ZBUF(I) ,1=1,12) 

C 

INOD2 = I NOD 
JNOD2 = JNOD 
KNOD2 = KNOD 
MARCH 2 = MARCH 
12 = 0 


END IF 


11 = 12 + 1 

12 = 12 + NSTORE 

READ! NU20 , 90 10 , ERR=400 ) ! XBUF ( I ) ,1=11,12) 
READ(NU20,9010,ERR=400) (YBUF( I) ,1=11,12) 
READ! NU20 , 9010 , ERR=400 ) ( ZBUF(I) ,1=11,12) 

GOTO 100 


200 CONTINUE 


WRITE (NU 6 ,9115) 

9115 FORMAT! ' END OF FILE REACHED ON FILE 20') 

NGRID = NGRID + 1 

IF( MARCH 2 . EQ. 1 )THEN 
ID! NGRID) = KNOD2 
JD( NGRID) = JNOD2 
KD( NGRID) = INOD2 
END IF 

IF! MARCH 2. EQ. 2 )THEN 
ID! NGRID) = INOD2 
JD( NGRID) = KNOD2 
KD( NGRID) = JNOD2 
END IF 

IF! MARCH 2 . EQ . 3 ) THEN 
ID! NGRID) = JNOD2 
JD( NGRID) = INOD2 
KD( NGRID) = KNOD2 
END IF 

ID! NGRID) = INOD2 
JD( NGRID) = JNOD2 
KD{ NGRID) = KNOD2 
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WRITE ( NU6 , 9 120 ) NGRID, ID{ NGRID) , JD< NGRID) , KD( NGRID) , 12 
WRITE ( NU6 , 8000 ) INOD2 , JNOD2 , KNOD2 
8000 FORMAT ( ' INOD2 , JNOD2 ,KNOD2 ',315) 

C 

WRITE ( 19,9010, ERR=450 ) 

1 ( XBUF (I), 1=1, 12), 

2 (YBUF( I) ,1=1,12) , 

3 ( Z8UF ( I ) , 1=1 , 12 ) 

C 

WRITE (NU6, 8 100 ) 

8100' FORMAT ( ' GRID WRITTEN TO 19 FILE') 

C 

REWIND(19) 

C 

IF( NGRID . GT • 1 ) THEN 
WRITE (NU6, 8 300) NGRID 
8300 FORMAT ( ' NGRID ’,15) 

WRITE (NU21,ERR=475) NGRID 
END IF 
C 

WRITE ( NU6 , 8350 ) ( ID(N) , JD(N) ,KD(N) ,N=1, NGRID) 

8350 FORMAT ( ' IDIM, JDIM,KDIM ’,315) 

C 

WRITE ( NU21, ERR=475 ) ( ID( N) , JD( N) , KD( N) ,N=1, NGRID) 

C 

C 

DO 300 N=l, NGRID 

12 = ID(N)*JD(N)*KD(N) 

C 

READ( 19,9010,ERR=450) 

1 (XBUF(I) ,1=1,12) , 

2 (YBUF(I) ,1=1,12) , 

3 (ZBUF(I) ,1=1,12) 

C 

WRITE ( NU2 1 , ERR=475 ) 

1 (XBUF(I) ,1=1,12) , 

2 ' ( YBUF (I), 1=1, 12), 

3 (ZBUF(I) ,1=1,12) 

C 

WRITE ( NU6 , 9130 ) N 

9130 FORMAT ( ' GRID ',15,' WRITTEN TO XYZHGM.DAT') 

C 

300 CONTINUE 
C 
C 

WRITE (NU6, 9 150) NGRID 

9150 FORMAT ( ' THERE ARE ',15,' GRIDS WRITTEN TO XYZHGM.DAT’) 
GOTO 500 
C 

400 WRITE ( NU6 , 9 140 ) 

9140 FORMAT( ' ERROR TRYING TO READ THE 20 FILE') 

GOTO 500 
C 

450 WRITE (NU6, 9 160) 

9160 FORMAT (' ERROR TRYING TO READ FILE 19’) 

GOTO 500 
C 

475 WRITE ( NU6,9170) 

9170 FORMAT ( ' ERROR TRYING TO WRITE' TO 20 FILE') 

C 
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C 


500 CONTINUE 
RETURN 


C 


END 


SUBROUTINE PROJECT (POINT, RANG E , ANG LE ) 


CALCULATES PROJECTION AND SWEEP ANGLE OF POINT ON BOWL SURFACE 
POINT = POSITION ON BOWL SURFACE 

RANGE = PROJECTION ONTO X AXIS FROM EDGE OF BOWL 
ANGLE = SWEEP ANGLE FROM VERTICAL 


DIMENSION POINT( 3 ) , BOWL( 3 ) , AXIS ( 3 ) ,VECTOR(3) ,SWEEP(3) ,VDUM(3) 

DATA BOWL /6 . 07 , 0 . 0 , 0 . 0/ 

DATA AXIS /l. 0,0. 0,0.0/ 

— DISTANCE ALONG AXIS FROM EDGE OF BOWL 

ONE = 1.0 

CALL VADD(ONE, POINT, -ONE, BOWL, VECTOR, VDUM) 

CALL VDOT( VECTOR, AXIS, RANGE) 

ANGLE FROM VERTICAL 

CALL VADD( ONE , VECTOR, -RANGE , AXIS , VDUM , SWEEP ) 

ANGLE = ACOS ( SWEEP ( 2 ) ) 

RETURN 

END 

*********************************************************************** 

* * * * * qt i LITY* ********************************************************** 
*********************************************************************** 

SUBROUTINE CROSS ( A , B ,C , N) 
n.C = CROSS PRODUCT OF A AND B (UNIT VECTOR) 


C 

C 


COMMON /COUNTER/ 
COMMON /HEADER/ 
COMMON /INITA/ 
COMMON /IWRITE/ 
COMMON /UNITS/ 
COMMON /ZONING/ 


NODES AV, NODETOT , NBNODES , NPLANE 

ITITLE ( 20 ) , LINE 

IDIM,MAPTEN, INCHES 

IWRTC , IWRTI , IWRTN 

NU5 , NU6 , NU19 , NU20 , NU2 1 

ISECT,NSECT, IZONE , NMBRNDS ( 3 ) 


DIMENSION A(3),B(3),C(3) 


Q-U7 


LOCKHEED-HUNTSVILLE ENGINEERING CENTER 



on n non o o n o o o o onoonoonoo o n non o non on 


LMSC-HEC TR F225755 


— CROSS PRODUCT 


C(l) = A( 2) *B{ 3 ) - A{ 3 ) *B( 2 ) 
C ( 2 ) = A( 3 ) *B ( 1 ) - A( 1 ) *B( 3 ) 
C( 3 ) = A( 1 ) *B ( 2 ) - A( 2 ) *B< 1 ) 

— MAGNITUDE 


CALL VMAG(C,CMAG) 

I F ( CMAG . GT .0.0) THEN 

NORMALIZE 

C(l) = C( 1 )/CMAG 
C( 2 ) = C(2)/CMAG 
C( 3 ) = C( 3 )/CMAG 


ELSE 

C(l) = 0.0 
C ( 2 ) = 0.0 
C( 3 ) = 0.0 

IF ( IWRTI . EQ . 0 ) RETURN 


ERROR PRINT (CROSS PRODUCT EQUALS ZERO) 


N * 10 **> 21 
30 -> 32 
40 -> 45 
50 -> 56 
60 -> 63 
70 -> 73 


SUBROUTINE EDGMAPl 
" EDGMAP2 

• SURFACE 

» BC 

* HOLES 

" FAIRING 


NODNUM = NODESAV + NODETOT 
I F ( N . LT .50 .OR. N.GT.59) NODNUM = NODNUM + 1 

LINE = LINE + 1 

IF( LINE .GE . 60 ) THEN 

WRITE(NU6, 1000) ITITLE, ISECT,NSECT, IZONE 
LINE = 2 

END IF 

WRITE(NU6, 1010) N, NODNUM 
END IF 


RETURN 

— FORMAT STATEMENTS 

000 FORMAT ( 1H1,10X,20A4,13X,8H SECTION , 12 , 3H OF,I3,9H FOR ZONE, 13) 

010 FORMAT ( 9H LOCATION , 1 3 , 36H s CROSS PRODUCT EQUALS ZERO FOR NODE, 16) 

END • 
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C 

SUBROUTINE ETABC ( MARCH , INDEX , NODE ) 


THIS ROUTINE CACULATES THE VALUE OF ETA. 


COMMON /INITA/ IDIM , MAPTEN , INCHES 

COMMON /INITC/ PI , RADDEG 

COMMON /INPUTA/ EDGE ( 6 , 1 2 ) , POINT ( 6 , 8 ) , S I DE ( 6 , 6 ) 

COMMON /INPUTBC/ ISIDE(3) 

COMMON /SPACING/ ISTRTCH(3) ,STRETCH(3) , ETAS ( 3 , 200 ) , ETA( 3 ) , DETA( 3 ) 
COMMON /ZONING/ ISECT , NSECT , I ZONE , NMBRNDS ( 3 ) 


FIRST NODE 


IF ( NODE . EQ . 1 ) THEN 

ETA( INDEX) = 0.0 

DETERMINE SIDE 

ISIDE ( INDEX) = 10 + (INDEX - 6)*INDEX 

STORE SPACING 

ETAS( INDEX, 1) = 0.0 
RETURN 
END IF 

LAST NODE 

IF ( NODE. EQ. NMBRNDS ( INDEX) ) THEN 

ETA( INDEX) = 1.0 

DETERMINE SIDE 

ISIDE{ INDEX) = 9 + (INDEX - 7)*INDEX/2 

— STORE SPACING 


ETAS (INDEX, NODE) = 1.0 

RETURN 

END IF 


INTERIOR NODES 


IF( ISECT. GT.l .AND. INDEX. NE .MARCH) GO TO 310 
— CALCULATE ETA 
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C 


ISTR = ISTRTCH( INDEX) + 1 

GO TO (100,110,120,130,140,150,160,170,180,190,200) ISTR 

EQUAL SPACING 

100 ETA( INDEX) = ETA( INDEX) + DETA( INDEX) 

GO TO 300 

INPUT ETA SPACING 

110 ETA( INDEX) = ETAS ( INDEX , NODE ) 

GO TO 300 

DECREASING SPACING; INPUT STRETCHING FACTOR 

120 RATIO = REAL( NODE - 1 ) /REAL( NMBRNDS ( INDEX ) - 1) 

XI = RATIO*STRETCH( INDEX ) /2 .0 
ETA1 = (EXP(Xl) - EXP(-Xl) )/( EXP( XI ) + EXP(-Xl)) 

X2 = STRETCH( INDEX)/2.0 

ETA2 = ( EXP ( X2 ) - EXP ( -X2 ) ) / ( EXP( X2 ) + EXP(-X2)) 

ETA( INDEX) = ETA1/ETA2 
GO TO 300 

INCREASING SPACING; INPUT STRETCHING FACTOR 

130 RATIO = REAL(NMBRNDS( INDEX) - NODE ) /REAL ( NMBRNDS ( INDEX) 

XI = RATIO* STRETCH (INDEX)/2.0 
ETA1 = ( EXP ( XI ) - EXP(-Xl) )/(EXP(Xl) + EXP(-Xl)) 

X2 = STRETCH( INDEX)/2.0 

ETA2 = (EXP(X2) - EXP( -X2 ) )/( EXP( X2 ) + EXP(-X2)) 

ETA( INDEX) = 1.0 - ETA1/ETA2 
GO TO 300 

DOUBLE STRETCHING; INPUT STRETCHING FACTOR — 

140 XI = ( STRETCH ( INDEX )/2 .0)/2 .0 

ETAMID = (EXP(Xl) - EXP ( -XI ) )/ ( EXP( XI ) + EXP(-Xl)) 

RATIO = REAL( NODE - 1 )/REAL( NMBRNDS ( INDEX) - 1) 

X2 = (RATIO - 0.5) * STRETCH ( INDEX )/2 . 0 
ETAMX = (EXP(X2) - EXP( -X2 ) ) / ( EXP( X2 ) + EXP(-X2)) 

ETA( INDEX) = (ETAMX + ETAMID)/ ( 2 . 0* ETAMI D) 

GO TO 300 

DECREASING SPACING; INPUT MINIMUM SPACING 


( 0 ) 


( 1 ) 


( 2 ) 


(3) 

1 ) 


(4) 


( 5 ) 
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150 ARGI = STRETCH( INDEX) *REAL( NODE - 1) 

EXPI = EXP(ARGI) 

EXP I I = 1.0/EXPI 

TANHI = (EXPI - EXPII)/ ( EXPI + EXPII) 

ARGN = STRETCH* INDEX) ‘REAL ( NMBRNDS ( INDEX) - 1) 

EXPN = EXP(ARGN) 

EXPNI = 1.0/EXPN 

TANHN = (EXPN - EXPNI)/(EXPN + EXPNI) 

ETA( INDEX) = TANHI/TANHN 
GO TO 300 

INCREASING SPACING; INPUT MINIMUM SPACING (6) 

160 ARGI = STRETCH* INDEX) *REAL( NMBRNDS* INDEX) - NODE) 

EXPI = EXP ( ARGI ) 

EXPII = 1.0/EXPI 

TANHI = (EXPI - EXPI I )/( EXPI + EXPII) 

ARGN = STRETCH ( INDEX) *REAL( NMBRNDS ( INDEX ) - 1) 

EXPN = EXP (ARGN) 

EXPNI = 1.0/EXPN 

TANHN = (EXPN - EXPNI )/( EXPN + EXPNI) 

ETA( INDEX) a 1.0 - TANHI/TANHN 
GO TO 300 

— DOUBLE STRETCHING; INPUT MINIMUM SPACING (7) 

170 ARGI = STRETCH( INDEX) *REAL( 2*NODE - NMBRNDS ( INDEX ) - 1) 

EXPI = EXP ( ARGI ) 

EXPII = 1.0/EXPI 

TANHI = (EXPI - EXPI I ) / ( EXPI + EXPII) 

ARGN = STRETCH( INDEX) ‘REAL ( NMBRNDS ( INDEX) - 1) 


EXPN = EXP (ARGN) 

EXPNI = 1.0/EXPN 

TANHN = (EXPN - EXPNI)/(EXPN + EXPNI) 

ETA( INDEX) = 0.5* (1.0 + TANHI/TANHN) 

GO TO 300 

— DECREASING SPACING; INPUT STRETCHING FACTOR (8) 

180 PIDN = PI/( STRETCH ( INDEX) ‘NMBRNDS ( INDEX) ) 

ETA( INDEX) = 1.0 - TAN ( PIDN* ( NMBRNDS ( INDEX) - NODE)) 

& /TAN( PIDN* (NMBRNDS{ INDEX) - 1)) 

GO TO 300 

— INCREASING SPACING; INPUT STRETCHING FACTOR (9) 

190 PIDN = PI/(STRETCH( INDEX) ‘NMBRNDS ( INDEX) ) 

C 

ETA( INDEX) = TAN { PIDN* ( NODE - 1)) 
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I 


& /TAN( PIDN* ( NMBRNDS ( INDEX) - 1)) 

GO TO 300 

USER INPUT STRETCHING FUNCTION (10) 

200 CONTINUE 

STORE SPACING 

300 ETAS ( INDEX, NODE) = ETA( INDEX) 

310 ETA( INDEX) = ETAS ( INDEX, NODE) 


RETURN 

END 

SUBROUTINE HOLES ( ANGLE , POINT , TANGENT, ICLOCK , IDUCT) 


CALCULATE POINT AND TANGENT ON HOLE 

ANGLE = ANGULAR LOCATION OF POINT 
POINT = POSITION AND FLOW VECTOR 
TANGENT = TANGENT 
ICLOCK = DIRECTION OF TANGENT 


COMMON /INITC/ PI , RADDEG 

DIMENSION POINT ( 6) ,TANGENT( 3) , US ( 3 ) , BC ( 3 ) , REF ( 3 ) 

ANG = ANGLE/RADDEG 

SANG = SIN(ANG) 

CANG = COS (ANG) 

IF( IDUCT. EQ. 1) THEN 

UPPER DUCT 


UNIT VECTOR ALONG AXIS 


A1 = 0.25506 

A3 = 0.092215 

A5 = 0.96252 


INTERMEDIATE VECTOR 


A2 

A4 

A6 


2 . 078 28*CANG - 
5.564 + 
-0.55072*CANG - 


0 . 05078 5*SANG 
2 . 1 4084 *SANG 
0 . 19165*SANG 


— REFERENCE POINT ON AXIS 

REF(l) = 9.445 
REF ( 2 ) = 5.564 
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REF ( 3 ) = 0.0 

ELSE 


LOWER DUCT 


UNIT VECTOR ALONG AXIS 

A1 = 0.40139 
A3 = 0.0 
A5 = 0.91591 

INTERMEDIATE VECTOR 

A2 = 1 . 96921*CANG - 2.11734 
A4 = 2 . 150*SANG 

A6 = 0.37334 - 0.86299*CANG 

REFERENCE POINT ON AXIS 

REF ( 1 ) = 7.32766 
REF( 2 ) = 0.0 
REF ( 3 ) = 0.37334 

END IF 


POSITION 


RAD = 7.890 

QUADRATIC EQUATION 

AA = A1*A1 + A3*A3 + A5*A5 

BB = 2 . 0* ( A1*A2 + A3*A4 + A5*A6) 

CC = A2*A2 + A4*A4 + A6*A6 - RAD**2 

ARG = BB*BB - 4.0*AA*CC 

— DISTANCE ALONG AXIS 

C - 0 . 5* ( -BB + SQRT(ARG) )/AA 

— VECTOR FROM SPHERE ORIGIN TO POINT ON HOLE 

US(1) » C*A1 + A2 
US( 2) = C*A3 + A4 
US ( 3 ) = C*A5 + A6 

— VECTOR FROM ORIGIN TO POINT ON HOLE 

POINT(l) = US(1) + 9.445 
POINT( 2 ) = US ( 2 ) 

POINT( 3 ) = US ( 3 ) 


TANGENT 


— VECTOR FROM AXIS TO POINT ON HOLE 

BC ( 1 ) = POINT(l) - C*A1 - REF ( 1 ) 
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BC( 2 ) = POINT( 2 ) - C* A3 - REF{ 2 ) 

BC( 3 ) = POINT( 3 ) - C*A5 - REF( 3 ) 

- — TANGENT AND FLOW VECTOR 

IF ( ICLOCK.EQ.O) THEN 

CALL CROSS(US,BC, TANGENT, 60) 

CALL CROSS { US , TANGENT, POINTf 4 ) ,61) 

ELSE 

CALL CROSS (BC, US, TANGENT, 62) 

CALL CROSS { TANGENT, US , POINT( 4 ) ,63) 

END IF 

RETURN 

END 

*********************************************************************** 

*****Q[j'fp|j'p**********************************************4************* 

A********************************************************************** 

SUBROUTINE OUTPUT( NUNIT, NODSTOR) 

PRINTOUT AND STORE DATA 

DOUBLE PRECISION NODE 
REAL NODE 

COMMON /COUNTER/ NODESAV, NODETOT, NBNODES , NPLANE 
COMMON /HEADER/ ITITLE ( 20 ) , LINE 
COMMON /INITA/ IDIM , MAPTEN , INCHES 
COMMON /IWRITE/ IWRTC , IWRTI , IWRTN 
COMMON /OUT/ NODE (5, 10000) 

COMMON /UNITS/ NU5 , NU6 , NU 19 , NU20 , NU21 
COMMON /ZONING/ ISECT , NSECT, IZONE, NMBRNDS ( 3 ) 

TOTAL NUMBER OF PLANES 

NPLANE = NPLANE + 1 

PRINT OUTPUT 

INTERMEDIATE PRINT 

100 IF( IWRTI. GT.O .AND. LINE. LE. 56) THEN 

WRITE(NU6, 1030) 

LINE = LINE + 3 
END IF 

— PRINT NODAL INFORMATION 
IPRINT = 1 
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C 

C 

C 

110 

C 

C 

C 

C 

C 


DO 120 1= 1 , NODSTOR 

IF( IWRTN.EQ.O .OR. I.LT.IPRINT) GO TO 120 

IPRINT = I + IWRTN 
NODNUM = I + NODETOT 

LINE = LINE + 1 

IF( LINE. GE. 60) THEN 

WRITE(NU6, 1000) ITITLE, ISECT,NSECT, I20NE 
WRITE(NU6, 1030) 

LINE = 5 

END IF 


WRITE(NU6, 1040) NODNUM, (NODE (J, I) ,J=1,5) 
120 CONTINUE 
STORE OUTPUT 


200 CALL BLKOUT(NUNIT, NODSTOR) 

—NODE COUNTERS 

NODETOT = NODSTOR + NODETOT 
NODESAV = NODESAV - NODSTOR 

— PRINT TOTAL NUMBER OF POINTS STORED 

IF(IWRTI.EQ.O) RETURN 

LINE = LINE + 1 


WRITE ( NU6 , 1060 ) NODSTOR, NPLANE , NUNIT , NODETOT 


RETURN 

FORMAT STATEMENTS 

1000 FORMAT ( 1H1 , 10X, 20A4, 13X, 8H SECTION, 12 , 3H OF,I3,9H FOR ZONE, 13) 

1020 FORMAT ( IX, 16 , 2 ( IX, F13 . 7 ) , 5X, F7 . 2 , 6X, F7 . 2 , 5X, 12 ) 

1030 FORMAT ( / 44H NODE X Y 2 ) 

1040 FORMAT { IX, 16 , 3 ( IX , F 1 3 . 7 ) , 3X, 2 ( 2X , F7 . 2 ) , 3X, 2 ( 2X , F7 . 2 ) , 5X, 1 2 ) 

1050 FORMAT ( 5X,I5,18H POINTS FROM PLANE, 14, 15H STORED ON UNIT, 13, 

& 23H: TOTAL POINTS STORED =,I6 ) 

1060 FORMAT ( 10X,I5,18H POINTS FROM PLANE, 14, 15H STORED ON UNIT, 13, 

S 23H: TOTAL POINTS STORED =,I6 ) 

C 

END 

C 

C*********************************************************************** 

q * * * * *OUTPUT ******* * **************************************************** 
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SUBROUTINE PICTURE( IDRAW) 


THIS ROUTINE DESCRIBES THE NOMENCLATURE 


COMMON /UNITS/ NU5 , NU6 , NU19 , NU20 , NU2 1 

WRITE ( NU6 / 300 ) 

WRITE ( NU6 ,310) 

WRITE(NU6,320) 

WRITE ( NU6 , 330 ) 

RETURN 

— FORMAT STATEMENTS 


300 FORMAT ( 

40X,37H 



3D NOMENCLATURE 



1 

/40X,19H 


IE 





2 

/40X,19H 


IT 





3 

/40X, 34H 


IA 

SURFACE 

4 



4 

/40X/ 32H 


12 

(TOP) 




5 

/40X, 50H 


POINT 4 I 




POINT 3 

6 

/40X,42H 


O- 



-O 


7 

/40X,42H 


8/1 

EDGE 3 

7/1 


8 

/40X, 50H 


E/ I 


E/ 

I 

S) 

310 FORMAT ( 

40X, 50H 


G/ I 

SURFACE 1 

G/ 

I 

U 

1 

/ 40X, 50H 


D/ I 

(BACK) 

D/ 

I 

E R 

2 

/40X, 50H 


E/ El 


E/ 

I 

D F 

3 

/40X, 50H 

POINT 8 O DI- 


— O POINT 

7 G A 

4 

/40X,50H 


I GI 

EDGE 11 

I 

I 

E C 

5 

/40X,50H 

S 

I El 


IE 

I 

E 

6 

/ 40X, 46H 

U 

El I 


ID 

I 

2 

7 

/ 40X, 50H 

R 

D I 41 

SURFACE 3 

IG 

I 

6 

8 

/40X,42H 

F 

G I I 

(FRONT) 

IE 

I 


9 

/40X, 50H 

A 

El O 


I 

0 

ETAl ) 

320 FORMAT ( 

40X, 51H 

C 

POINT 1 /5 

EDGE 1 

IT 

/ 

POINT 2 

1 

/40X, 42H 

E 

11 /E 


10 

/6 


2 

/40X, 41H 


2 I /G 


I /E 


3 

/40X, 40H 

5 

I /D 


I /G 



4 

/40X, 39H 


I /E 


I /D 



5 

/40X, 38H 


0/ 


O/E 



6 

/40X, 44H 

POINT 5 

EDGE 9 

POINT 

6 

7 

/40X,10H 


/ 





8 

/40X,26H 


/3 SURFACE 2 




9 

/40X, 25H 


/A (BOTTOM) 



) 

330 FORMAT ( 

40X, 8H 


/T 





& 

/40X, 7H 


/E 




) 


END 

SUBROUTINE VADD( CA, A, CB , B , C , UC ) 
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VADD COMPUTES C,THE SUM OF VECTORS CA*A AND CB*B, WHERE CA AND 
CB ARE SCALARS. UC IS A UNIT VECTOR DIRECTED ALONG C. 


DIMENSION A(3),8<3),C(3),UC(3) 
SUM = 0.0 
DO 10 1=1,3 

C(I) = CA*A( I ) + CB*B(I) 

10 SUM = SUM + C(I)*C(I) 

CM AG = SORT (SUM) 

RMAG = 0.0 

IF(CMAG.GT.O.O) RMAG = 1.0/CMAG 

UC(1) = C(1)*RMAG 
UC( 2) - C( 2 ) * RMAG 
UC( 3 ) = C ( 3 ) * RMAG 

RETURN 

END 

SUBROUTINE VDOT(A,B,C) 


VDOT COMPUTES C ,THE DOT PRODUCT OF VECTORS A AND B. 


DIMENSION A( 3 ) ,B ( 3 ) 

C = 0.0 

DO 10 1=1,3 
10 C = C + A ( I ) *B ( I ) 

RETURN 
END 

SUBROUTINE VMAG( VECTOR, VECMAG) 

VMAG DETERMINES THE MAGNITUDE OF A VECTOR 


DIMENSION VECTOR( 3 ) 

C 

VECMAG = SORT ( VECTOR( 1 ) * * 2 + VECTOR<2)**2 + VECTOR( 3)**2) 
C 

IF(VECMAG.LT. 0.0000001) VECMAG = 0.0 
C 

RETURN 
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END 

SUBROUTINE FAIRING( PT1 , PT2 , POINT, TANGENT, SNORMAL) 


TRANSITION FROM BOWL TO DUCT 

PT1 = POINT ON HOLE (ON SPHERE) 

PT2 = POINT ON DUCT (END OF FAIRING) 
POINT = POSITION 
TANGENT = TANGENT 
SNORMAL = SURFACE NORMAL 


COMMON /COUNTER/ NODESAV, NODETOT, NBNODES , NPLANE 
COMMON /INITC/ PI / RADDEG 

COMMON /SPACING/ ISTRTCH( 3) , STRETCH ( 3 ) , ETAS ( 3 , 200 ) , ETA( 3 ) , DETA( 3 ) 
COMMON /20NING/ ISECT, NSECT, I ZONE ,NMBRNDS ( 3 ) 

DIMENSION PT1 ( 3 ) , PT2 ( 3 ) , POINT( 6 ) ,TANGENT( 3 ) ,SNORMAL(3) 

DIMENSION PTW ( 3 ) ,UC(3) , CP( 3) ,REF( 3) , VDUM( 3) 

DIMENSION Rl(3) ,R2(3) ,UN(3) ,UP(3) 

NODNUM = NODESAV + NODETOT + 1 

—VECTOR FROM SPHERE ORIGIN TO EDGE OF HOLE 

PI = PT1(1) - 9.445 
P2 = PT 1 ( 2 ) 

P3 = PT 1 ( 3 ) 

PMAG = SQRT(P1*P1 + P2*P2 + P3*P3) 

— NORMALIZED 

SI = Pl/PMAG 
S 2 = P2/PMAG 
S3 = P3/PMAG 

IDUCT = IZONE-3 

WRITE ( 6 > 3030 ) IAXIS , JAXIS , KAXIS 
3030 FORMAT ( ' IAXIS( 2) ,JAXIS( 3) ,KAXIS( 1) =’,3110) 

I F ( IDUCT. EQ. 1) THEN 


UPPER DUCT 


—AXIS 


U1 = .0.25506 
U 2 = 0.092215 
U3 = 0.96252 

— INTERMEDIATE VALUES 

CADD = 5 . 564*U2 
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C 

A1 = Sl*(1.0 - 01**2) - S2*U1*U2 - S3*U1*U3 

A2 = S 2 * ( 1.0 - U 2 * * 2 ) - S1*U1*U2 - S3*U2*U3 

A3 = S3* ( 1 . 0 - U3**2) - S1*U1*03 - S2*U2*U3 

C 

B 1 = 5. 564*U1*U2 

B2 = 5 . 564* ( U2*U2 - 1.0) 

B3 = 5 . 564*U2*U3 

REFERENCE POINT ON AXIS 

REF(l) = 9.445 
REF( 2) = 5.564 
RE F ( 3 ) = 0.0 

WRITE ( 6 » 30 10 ) 

3010 FORMAT ( ' UPPER DUCT') 


ELSE 


LOWER DUCT 


AXIS 

U1 = 0.40139 
U2 = 0.0 
U3 = 0.91591 

— INTERMEDIATE VALUES 

CADD = -2 . 1 1734*U1 + 0.37334*U3 

A1 = Sl*(1.0 - U1*U1) - S3*Ul*U3 
A2 = S2 

A3 = S3* ( 1 . 0 - U3*U3 ) - S1*U1*U3 

B 1 = -2 . 1 1734* ( U1*U1 - 1.0) + 0. 37334*U1*U3 
B2 = 0.0 

B3 = -2. 11734*U1*U3 + 0 . 37334* { U3*U3 - 1.0) 

REFERENCE POINT ON AXIS 

REF(l) = 7.32766 
REF ( 2 ) = 0.0 
REF ( 3 ) = 0.37334 

WRITE( 6, 3020) 

3020 FORMAT ( ' LOWER DUCT’) 

END IF 


C 

C 


CALCULATE END POINT OF FAIRING 
—QUADRATIC EQUATION 

AA = A1*A1 + A2*A2 + A3*A3 - 1.0 

BB = 2 . 0* ( A1*B1 + A2*B2 + A3*B3 + 5.930) 

CC = B1*B1 + B2*B2 + B3*B3 - 35.1649 
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C 


DD = BB*B8 - 4 . 0* AA*CC 


C 

SDD = SQRT(DD) 

C 

RH01 = 0 . 5* ( -BB + SDD)/AA 
RH02 = 0 . 5* ( -BB - SDD)/AA 


RADIUS FROM SPHERE ORIGIN TO ARC CENTER 

RHO = RHOl 
I F ( RH02 . GT . RHOl ) RHO = RH02 

ARC RADIUS 


RAD = RHO - 7.890 

— RADIUS FROM DUCT CENTERLINE TO ARC CENTER 
CN = 1.96 + RAD 

— UNIT VECTOR FROM DUCT CENTERLINE TO ARC CENTER 

Cl = (RHO*Al + B1)/CN 
C2 a (RHO*A2 + B2)/CN 
C3 = ( RHO* A3 + B3)/CN 

SWEEP ANGLE OF ARC 

SN = C1*S1 + C2*S2 + C3*S3 

ANG = ACOS(SN) 

VECTOR FROM ORIGIN TO ARC CENTER 

DC1 = RHO*S 1 + 9.445 
DC 2 = RHO*S2 
DC 3 = RHO*S3 


— DISTANCE ALONG DUCT CENTERLINE 

CA = RHO* ( S1*U1 + S2*U2 + S3*U3) - CADD 
—VECTOR FROM ORIGIN TO END POINT OF FAIRING 


PTW(l) = CA*U1 + 1.96*C1 + REF(l) 
PTW( 2 ) = CA*U2 + 1 . 96*C2 + REF ( 2) 
PTW( 3 ) = CA*U3 + 1 . 96*C3 + REF( 3) 


WRITE ( 6 , 3000 ) 
3000 FORMAT (' NODE 
1 
2 


ANGL = ANG*RADDEG 

NODNUM,RHO,RAD,CA,ANGL,Cl,C2,C3, (PTW( I) ,1=1,3) 
=',110 / ' RHO, RAD, CA, ANGL = ' , 4 ( 2X, F12 . 5 ) 

/ ' Cl ,C2,C3 =' , 3 ( 2X , F12 . 5 ) 

/ ' PTW =' , 3 ( 2X, F12 . 5 ) ) 


JUNCTION = 7 


IF(IAXIS.LT. JUNCTION) THEN 
FAIRING SURFACE 
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VECTOR FROM ARC CENTER TO HOLE 

Rl(l) = PT1(1) - DC1 
R1 ( 2 ) = PT1 ( 2 ) - DC2 
R1 ( 3 ) = PT1 ( 3 ) - DC 3 

VECTOR FROM ARC CENTER TO END OF FAIRING 

R2 ( 1 ) = PTW(l) - DC1 
R2 ( 2 ) = PTW( 2 ) - DC2 
R2 ( 3 ) = PTW( 3 ) - DC 3 

LOCAL COORDINATE SYSTEM 

CALL CROSS (R1,R2,UN,70) 

CALL CROSS ( RI / UN , UP/ 7 1 ) 

CALL CR0SS(UN,UP,R1,72) 

CALL VDOT(UP,R2,UPE) 

IF(UPE.LT.O.O) CALL CROSS ( UN , Rl , UP , 73 ) 

ARC ANGLE ALONG FAIRING 

THET = ANG* ( FLOAT ( IAXIS ) - 1.0)/3.0 

CANG = COS (THET) 

SANG = SIN(THET) 

UNIT VECTOR FROM ARC CENTER TO POINT ON FAIRING 

CALL VADD( CANG , Rl , SANG , UP , VDUM , UC ) 

POSITION 

POINT ( 1 ) = DC 1 + RAD*UC( 1 ) 

POINT ( 2 ) « DC 2 + RAD*UC( 2) 

POINT ( 3) = DC 3 + RAD*UC( 3) 

WRITE (6,3050) PT1 , PT2 , RHOl , RH02 , POINT 
3050 FORMAT ( ' PT1,PT2 = ' , 6 ( 2X, F12 . 5 ) 

1 / ' RHOl, RH02, POINT =' , I 4X, 5 ( 2X, F12 . 5)) 

TANGENT 

CALL VADD(CANG, UP, -SANG, Rl , VDUM, TANGENT) 

C 

CE = TANGENT( 1 ) *U1 + TANGENT( 2 ) *U2 + TANGENT( 3 ) *U3 
C 

IF(CE.LT.O.O) THEN 
C 

TANGENT ( 1 ) => -TANGENT ( 1 ) 

TANGENT ( 2 ) = -TANGENT { 2 ) 

TANGENT ( 3 ) = -TANGENT ( 3 ) 
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— NORMAL 

SNORMAL(l) = — UC ( 1 ) 
SNORMAL( 2 ) = -UC ( 2) 
SNORMAL( 3 ) = -UC(3) 

ELSE 


DUCT SURFACE 


RATIO BETWEEN END OF FAIRING AND END OF SECTION 

DUCT SURFACE 


RATIO BETWEEN END OF FAIRING AND END OF SECTION 

TOTND = FLOAT( NMBRNDS ( 2) ) 

RATIO = (TOTND - IAXIS) / (TOTND - JUNCION) 

STR21 = 3.0 

XI = RATIO* STR2 1/2.0 

ETA1 = (EXP(Xl) - EXP(-Xl) )/( EXP(X1 ) + EXP(-Xl)) 
X2 = STR21/2.0 

ETA2 = (EXP(X2) - EXP(-X2) )/(EXP(X2) + EXP(-X2)) 
EPS = 1.0 - ETA1/ETA2 
EPS 1 = 1.0 - EPS 

VECTOR FROM ORIGIN TO POINT ON DUCT 

CALL VADD( EPS 1 , PTW, EPS , PT2 , CP , VDUM ) 

POSITION 

POINT ( 1 ) = CP(1) 

POINT( 2) = CP( 2 ) 

POINT( 3 ) = CP( 3 ) 

TANGENT 

TANGENT(l) » U1 
TANGENT ( 2) = U2 
TANGENT ( 3 ) = U3 

NORMAL 

SNORMAL(l) = Cl 
SNORMAL( 2) = C2 
SNORMAL( 3 ) = C3 

END IF 


FLOW VECTOR 


POINT( 4 ) = TANGENT ( 1 ) 
POINT( 5 ) = TANGENT ( 2) 
POINT( 6) = TANGENT ( 3 ) 
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C 

RETURN 

END 

C 

C**EOR** 

C**EOI** 
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